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1. Introduction 

In the past two decades, next-to-leading order (NLO) QCD computations have become 
standard tools for phenomenological studies at lepton and hadron colliders. QCD tests 
have been mainly performed by comparing NLO results with experimental measurements, 
with the latter corrected for detector effects. 

On the experimental side, leading order (LO) calculations, implemented in the context 
of general purpose Shower Monte Carlo (SMC) programs, have been the main tools used 
in the analysis. SMC programs include dominant QCD effects at the leading logarithmic 
level, but do not enforce NLO accuracy. These programs were routinely used to simulate 
background processes and signals in physics searches. When a precision measurement was 
needed, to be compared with an NLO calculation, one could not directly compare the 
experimental output with the SMC output, since the SMC does not have the required 
accuracy. The SMC output was used in this case to correct the measurement for detector 
effects, and the corrected result was compared to the NLO calculation. 

In view of the positive experience with QCD tests at the NLO level, it has become 
clear that SMC programs should be improved, when possible, with NLO results. In this 
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way a large amount of the acquired knowledge on QCD corrections would be made directly 
available to the experimentalists in a flexible form that they could easily use for simulations. 

The problem of merging NLO calculations with parton shower simulations is basically 
that of avoiding overcounting, since the SMC programs do implement approximate NLO 
corrections already. Several proposals have appeared in the literature [1, 2, 3, 4] that can be 
applied to both e + e~ and hadronic collisions, and two approaches [5, 6] suitable for e + e~ 
annihilation. Furthermore, proposals for new shower algorithms, that should be better 
suited for merging with NLO results, have appeared in the literature (see refs. [7, 8, 9, 10, 
11, 12]). 

The MCONLO proposal [2] was the first one to give an acceptable solution to the over- 
counting problem. The generality of the method has been explicitly proven by its applica- 
tion to processes of increasing complexity, such as heavy-flavour-pair [13] and single-top [14] 
production. 1 The basic idea of MC@NL0 is that of avoiding the overcounting by subtracting 
from the exact NLO cross section its approximation, as implemented in the SMC pro- 
gram to which the NLO computation is matched. Such approximated cross section (which 
is the sum of what have been denoted in [2] as MC subtraction terms) is computed an- 
alytically, and is SMC dependent. On the other hand, the MC subtraction terms are 
process-independent, and thus, for a given SMC, can be computed once and for all. In 
the current version of the MCONLO code, the MC subtraction terms have been computed 
for HERWIG [15]. In general, the exact NLO cross section minus the MC subtraction terms 
does not need to be positive. Therefore MC@NL0 can generate events with negative weights. 
For the processes implemented so far, negative-weighted events are about 10-15% of the 
total. Their presence does not imply a negative cross section, since at the end physical 
distributions must turn out to be positive. 

The features implemented in MC0NL0 can be summarized as follows: 

- Infrared-safe observables have NLO accuracy. 

- Collinear emissions are summed at the leading-logarithmic level. 

- The double logarithmic region (i.e. soft and collinear gluon emission) is treated cor- 
rectly if the SMC code used for showering has this capability. 

In the case of HERWIG this last requirement is satisfied, owing to the fact that its shower is 
based upon an angular-ordered branching. 

In ref. [4] a method, to be called POWHEG in the following (for Positive Weight Hard- 
est Emission Generator), was proposed that overcomes the problem of negative weighted 
events, and that is not SMC specific. In the POWHEG method the hardest radiation is gen- 
erated first, with a technique that yields only positive- weighted events using the exact NLO 
matrix elements. The POWHEG output can then be interfaced to any SMC program that 
is either p T -ordered, or allows the implementation of a p T veto. 2 However, when interfacing 

X A complete list of processes implemented in MC@NL0 can be found at 

http : / /www . hep . phy . cam . ac . uk/ theory/webber/MCatNLO . 

2 A11 SMC programs compatible with the Les Houches Interface for User Processes [16] should comply 
with this requirement. 
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POWHEG to angular-ordered SMC programs, the double-log accuracy of the SMC is not 
sufficient to guarantee the double- log accuracy of the whole result. Some extra soft radia- 
tion (technically called vetoed-truncated shower in ref. [4]) must also be included in order 
to recover double-log accuracy. In fact, angular ordered SMC programs may generate soft 
radiation before generating the radiation with the largest p T , while POWHEG generates 
it first. When POWHEG is interfaced to shower programs that use transverse- momentum 
ordering, the double-log accuracy should be correctly retained if the SMC is double-log ac- 
curate. The ARIADNE program [17] and PYTHIA 6.4 [18] (when used with the new showering 
formalism), both adopt transverse- momentum ordering, in the framework of dipole-shower 
algorithm [19, 20, 21], and aim to have accurate soft resummation approaches, at least in 
the large N c limit (where N c is the number of colours). 

A proof of concept for the POWHEG method has been given in ref. [22], for ZZ produc- 
tion in hadronic collisions. In ref. [23] the method was also applied to QQ hadroproduction. 
Detailed comparisons have been carried out between the POWHEG and MC0NL0 results, 
and reasonable agreement has been found, which nicely confirms the validity of both ap- 
proaches. In ref. [24] the POWHEG method, interfaced to the HERWIG Monte Carlo, has 
been applied to e + e" annihilation, and compared to LEP data. The method yields better 
fits compared to HERWIG with matrix-element corrections. The authors of ref. [24] have also 
provided an estimate of the effects of the truncated shower, which turned out to be small. 

In the present work we give a detailed description of the POWHEG method. Our 
aim is to provide all of the necessary formulae and procedures for its application to general 
NLO calculations. We first formulate POWHEG in a general subtraction scheme. Then, we 
illustrate it in detail in two such schemes: the Frixione, Kunszt and Signer (FKS) [25, 26] 
and the Catani and Seymour (CS) [27] one. The CS method has been widely used in 
the literature. On the other hand, the FKS method has already been used extensively in 
the MC0NL0 implementations. Furthermore, the NLO cross sections for vector-boson and 
heavy-quark pair production used in the POWHEG implementations of refs. [22, 23] have 
a treatment of initial-state radiation that is essentially the same one used in FKS. 

Our paper is organized as follows. In section 2 we summarize the general features of 
the NLO computations and of the subtraction formalisms. In section 2.4 all the details 
of the FKS subtraction method are given, and in section 2.5 the basic features of the CS 
approach are summarized. 

In section 3 a general discussion of the inclusion of NLO corrections in a parton shower 
framework is given, together with a basic introduction of the POWHEG method. 

In section 4 we go through all the details of the POWHEG method. The method is 
presented in general, and it is shown how to apply it within any subtraction framework. 
Thus, this section does not refer in particular to either the FKS or the CS method. 

In section 4.4 we discuss the accuracy of the POWHEG approach in the resummation of 
soft-gluons effects. We show that, with an appropriate prescription for the evaluation of the 
running coupling used in POWHEG for the generation of radiation, one can easily obtain 
next-to-leading logarithmic (NLL) accuracy in soft-gluon radiation, provided the process 
in question has no more than three incoming or outgoing coloured partons at the Born 
level. If the Born process involves more than three coloured partons, there are left-over 
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soft terms that are not correctly represented by POWHEG. We show, however, that with 
a modest modification of the algorithm one can also correctly resum these contributions at 
the level of the leading terms in a large N c expansion. 

In section 5 the formulae needed for the explicit construction of a POWHEG in the 
FKS framework is given. The same is done in section 6 for the CS method. 

Finally, in section 7 two simple examples are discussed in both the FKS and the CS 
framework, namely the production of hadrons in e + e" annihilation, and the production of 
a massive vector boson (or a virtual photon) in hadronic collisions. 

This paper is considerably long, and it involves many technical details. The length is 
partly due to the fact that we deal with two subtraction methods. The reader may skip the 
one she/he is not interested in. The example sections are particularly long and pedantic. 
The reader may not be interested in reading all of them. Section 4.4 is technically complex, 
but it may be almost completely skipped on a first reading. 

2. NLO computations 
2.1 Generalities 

In this section, we describe the general features of an NLO calculation for a generic hadron- 
hadron collision process. In lepton-hadron and lepton-lepton collisions, the treatment is 
similar, but simpler. For example, in the case of lepton-lepton collisions, the parton- 
distribution functions for the incoming particles are replaced by delta functions. 

We consider 2 — > n processes, where the momenta of the particles satisfy the momen- 
tum conservation 

x (B K m + x e K e = k 1 + ... + k n , (2.1) 

where x® are the momentum fractions of the incoming partons, and K @ the four-momenta 
of the incoming hadrons. In what follows, we also use the notation 

to denote the momenta of the incoming partons. We define, as usual, 

S = (K 9 + K e f , s = (k 9 + k e f . (2.3) 

We denote by 3>„ the set of variables 

*„ = {x ffi ,x e , fei, . . . , k n } (2.4) 

constrained by momentum conservation (eq. (2.1)), and by the on-shell conditions for final- 
state particles. We collectively denote by B the squared matrix elements 3 relevant to the 
LO contributions to our process. The total cross section at leading order is given by 

a hQ = J d® n CB{® n ) (2.5) 

3 We always assume spin and colour sums and averages when needed, and the inclusion of the appropriate 
flux factor. 
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where C is the parton luminosity 4 



C = C(x (B ,x e ) = / e (x e )/ e (x e ), (2.6) 

and 

d& n = dx 9 dx e d$ n (fe e + /c e ; fei, . . . , k n ) , (2.7) 
with d$ n the n-body phase space 

/ n \ n , 3 , 

d<t> n (q; h, . . . , fc n ) = (2vr) 4 <5 4 U - £ fcij [J • ( 2 - 8 ) 

In case of leptons in the initial state, the corresponding parton distribution function f{x) 
in eq. (2.6) is replaced by 5(1 — x). 

The real contributions at the NLO arise from the tree-level squared amplitudes for 
the 2 — > n + 1 parton process, which we denote by 1Z. As before, we denote by 3? n +i the 
corresponding set of variables 

*n+i = {x s ,x e ,ki, . . . ,k n+1 } (2.9) 

constrained by momentum conservation and on-shell conditions. 

The virtual contributions arise from the interference of the one-loop amplitudes times 
the LO amplitudes. We denote by Vb the renormalized virtual corrections, that is, we 
assume that all ultraviolet divergences have already been removed by renormalization. 
These terms still contain infrared divergences. Therefore, they are computed in d = 4 — 2e 
dimensions, and the divergences appear as 1/e 2 and 1/e poles. The subscript b (for "bare") 
reminds us of the presence of infrared divergences in the amplitude. 

In hadronic collisions, the complete cancellation of the initial-state collinear singular- 
ities is achieved by adding two counterterms, one for each of the incoming partons (©, 
©), to the differential cross section. We denote them by and G e ,b- The factorization 
counterterms are infrared divergent in four dimensions. Therefore, they are computed in 
d = 4 — 2e dimensions, and the divergences appear as 1/e poles. To remind this fact, also 
in this subscript b has been included in the notation. 

The total NLO cross section is given by 5 

^nlo = J d® n C [£($„) + V b (* n )] + J d® n+1 £K(® n+1 ) 

+ J d* n , ffi £S ffiib (* n , ffl ) + J d$ n , e £S eib ($ n , e ) , (2.10) 



where 



d* n +i = dx 9 dx e d$ n+1 (fc e + fc e ; fei, . . . , k n+1 ) . (2.11) 



4 In this section we drop the parton flavours and the scale dependence in the luminosity, for ease of 
notation. 

5 The (5©,b terms are present only for incoming hadrons. If one or both the incoming particles are leptons, 
the corresponding Qb is zero. 
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^n,@ denotes configurations in which one of the final-state partons is collinear to one of the 
incoming partons. Thus, such configurations are effectively n-body final-state ones, except 
for the energy degree of freedom of the parton collinear to the beam. We then write 



^n,e — { x ®i x ei z i • • • > ^n} ) z x s>K® + x qKq — ^ ] ki i (2.12) 

i=\ 
n 

^n,e = { x ®i x qi z i ki, ■ ■ ■ , k n } , x®K 9 + zXqKq = ^ ^ ki , (2.13) 



i=i 



where z is the fraction of momentum of the incoming parton after radiation. We can 
associate with the phase-space configuration <& n ,® 

an underlying n-body configuration *l* n 

defined as 

&n = { x ®, x ei k\, . . . k n } , x© = zx & , 5g = x^ . (2-14) 

Thus, the values of x© in the underlying n-body configuration are constrained by momen- 
tum conservation, and do not depend upon z. We also define 

d&n,® = dx @ dx e dz d$ n (z k @ + k e ; h, . . . , k n ) , (2.15) 
d<& n ,e = dx 9 dx e dz d& n (k @ + z k e ; ki, . . . , k n ) . ( 2 -16) 

We now consider a generic observable O, function of the final-state momenta. O could 
be, for example, a product of theta functions describing a particular histogram bin for the 
distribution of some kinematic observable. Its expected value is given by 

(O) = J d& n £O n (3> n ) [#(<!>„) + V b (<I> n )" 

+ J d$ n+1 £O n+1 ($ n+1 ) ft(* n+ i) 

+ J d*„, e £O n (* n ) Se,b(*n,e)+ / d$ n , e £O n (® n ) a e ,b(*n,e), (2-17) 

where O n and O n +\ are the expressions of the observable O in terms of n and (n + 1) final- 
state particle momenta, and 3> n , in the <&«,© integrals, is the corresponding underlying 
n-body configuration. We require that O is an infrared-safe observable, and, furthermore, 
we require that the Born contribution in eq. (2.17) (i.e. the term proportional to B) is 
infrared finite (thus, for example, if our n-body process corresponds to Z + jet production, 
the observable O n must suppress the region where the jet is emitted at low transverse mo- 
mentum). Under these assumptions, the real matrix elements contribution (i.e. the term 
proportional to TV) is finite in the whole phase space d<& n+ i, except for the regions that 
correspond to soft and collinear emissions. There, the divergences are integrable only in d 
dimensions, and yield 1/e 2 and 1/e poles. Furthermore the divergences of each term on the 
r.h.s. of eq. (2.17) cancel in the sum, and the total cross section is finite. Observe that the 
argument of O in the last two terms on the right hand side of eq. (2.17) is set equal to $ n 
rather than 3? ra ,©> owing to the fact that O is an infrared-safe observable. The integrals 
in eq. (2.17) are usually too difficult to be performed analytically (because of the involved 
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functional form of O) and, being divergent, they are not suited for numerical computations. 
For these reasons, different strategies have been proposed for the computation of observ- 
ables in QCD. One of the most successful is the so-called subtraction method, pioneered 
in ref. [28], which we discuss in the next section. 

2.2 Subtraction formalism 

The subtraction formalism requires the definition of a set of functions C^ a \ called real 
counterterms. Each a is associated with a particular singular region, i.e. with a configura- 
tion that has either a final-state parton with four momentum equal to zero, or a final-state 
massless parton with momentum proportional to an initial-state or to another final-state 
massless parton. Furthermore, for each a, a mapping 6 M*") 

?i = M< a »(M , *& = {&\&\kf\...,~k$ 1 } (2-18) 

is defined that maps the (n + l)-body configuration into a singular one. 

The real counterterms and the mapping have the following property: for any infrared- 
safe observable O n+ i(<& n+ i), that vanishes fast enough if 3> n +i approaches two singular 
regions at the same time, the function 

ft(*„+i) 0„+i(*n+i) - £c< a >(* n+1 ) O n+1 (mW (* n+1 )) (2.19) 

a 

has at most integrable singularities in the <fr n +i space. Observe that the above condition 
does not always imply that 

ft(* n+1 )-^C^(*n+i) (2.20) 

a 

is also integrable. This is the case if the corresponding n-body process has no singularities, 
like, for example, in Z production in hadronic collisions. 

Each singular region a is characterized by a different mapping, and, for this reason, 
we use the superscript a on the tilded variables. For ease of notation, we use the following 
context convention: if an expression is enclosed in the subscripted squared brackets 

[-]«> ( 2 - 21 ) 

we mean that all variables appearing inside have, when applicable, the superscripts corre- 
sponding to the subscript of the bracket. Thus we write 



n+l 



|x e ,£ e ,fei,...,fe n+ ijj . (2.22) 



The form of the singular configurations ^^+1 differs according to the nature of the singular 
region. More specifically: 



6 In some approaches, the counterterms are different from zero only in a finite neighborhood of the 
corresponding singular regions. In these cases, the mapping needs to be defined only there. 
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• If a is associated with a soft (S) region, the singular configuration has a final-state 
parton with null four-momentum. 

• If a is associated with a final-state collinear singularity (FSC), the singular configu- 
ration has two massless final-state partons with parallel three-momenta. 

• If a is associated with an initial-state collinear (ISC) singularity, the singular con- 
figuration has a massless outgoing parton with three-momentum parallel to the mo- 
mentum of one incoming parton. 

The mapping (2.18) must be smooth near the singular region, and it must become the 
identity there. In other words if, for instance, a is associated with the FSC region where 
the particles i and j become collinear, we must have ^"+1 = ^n+i for &i II kj. Notice that 
also the x's of the configurations ^^+1 do not necessarily coincide with x ffi and x e for all 
(n + l)-body configurations <& n +i- On the other hand, they do coincide in the singular 
limit. 

As in the case of the collinear configurations (eqs. (2.12) and (2.13)), we associate with 
each 3^+! configuration an n-body configuration that we will call the underlying 

n-body configuration 

= [{x 9 ,x e ,h,...,k n }] a • (2-23) 

is obtained as follows: 

• If a G S (i.e. it is a soft region), <j>£^ is obtained by deleting the zero momentum 
parton. 

• If a G FSC (i.e. it is a final-state collinear region), <j>j^ is obtained by replacing the 
momenta of the two collinear partons with their sum. 

• If a G ISC (i.e. it is an initial state collinear region), is obtained by deleting the 
radiated collinear parton, and by replacing the momentum fraction of the initial-state 
radiating parton with its momentum fraction after radiation. 

In all the above cases, the final-state momenta are relabelled with an index that takes 
values in the range 1, . . . ,n. Observe that, as a consequence of the procedure itemized 
above, the variables in <j>£^ are constrained by momentum conservation 

n 

x e K (B +x e K e = J2k ■ ( 2 - 24 ) 
Furthermore, for S or FSC regions, we have 

= £© • (2.25) 
This does not hold for ISC regions: in the © direction, for example, we have 

^© < > = > (2.26) 
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and the analogous one for the case of ISC in the © direction. 

We stress that the difference in our notation between the and is a minor 

one: the former has an unresolved parton, while in the latter all partons are resolved. On 
the other hand, it is necessary to introduce (together with the concept of underlying 
n-body configuration), since it is formally the argument of Born-like matrix elements, and 
will play a central role in the development of the POWHEG formalism. 

In the subtraction method one rewrites the contribution to any observable O coming 
from real radiation in the following way 

J d3> n+1 £O n+1 ($ n+1 ) K{® n+l ) = Y J j d® n +i [C0 n (& n ) C(* n+ i)] a + 

J d$ n+l {cO n+l {3> n+l )1l{3> n+l )-Y J [£O n {3> n ) C(* n+ i)] Q }, (2.27) 

a 

where C = £(x e ,x e ). In this way, under the assumptions we have made about the coun- 
terterms, and the assumption that O is an infrared-safe observable, the second term on the 
r.h.s. of eq. (2.27) is integrable in d = 4 dimensions. 

The first term on the r.h.s. of eq. (2.27) is divergent. In order to deal with it, we 
introduce, for each a, the (n + l)-phase space parametrization 



71+1 



and the corresponding phase-space element 



(2.28) 



(2.29) 



In words, we parametrize the (n + l)-phase space in terms of an n-body phase space 
(obtained as described earlier), plus (three) more variables that describe the radiation 
process. The left-right arrow in eq. (2.28) indicates that the correspondence is one to one. 7 
The range of the radiation variables in ^^ad 

may depend upon $1"'. Furthermore, eq. (2.29) 
implicitly defines a Jacobian, possibly dependent upon «&^*\ that we conventionally include 
into d^rad- We- call the parametrization (2.28) the emission factorization. 

We now distinguish two cases: the FSC+S case and the ISC one. In the former case 
we have 

t = £(x e , x e ) = £(x ffi , x e ) . (2.30) 



Defining 

C(* n ) = f d$ rad C(* n+ i) 

J -I aG{FSC,s} 

we can write the generic term in the first sum on the r.h.s. of eq. (2.27) as follows 



d& n+1 £O n (& n ) C(* n+ i) = / d^ n CO n (^ n )C{^ n ) 



- ae{Fsc,s} 



(2.31) 



(2.32) 



7 The correspondence needs only to be defined where the corresponding counterterm is non- vanishing. 
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In the ISC case, we cannot factor out the luminosity so easily, since x @ / x@. We define 



C(*„,z) = y d$ rad C($ n+ i) zS(z-x @ / 

which formally introduces the momentum fraction z, and write 

j d^ n+1 £O n (^ n ) C(* n+ i) = y d* n ^£O n (* n )C (*„,*) 

Notice that, owing to the delta function in eq. (2.33) we have 

r - r(~ ~ \ - / £( s ffi/ z '^e) for a G ISC a 
I £(x®,x e /z) for q G ISC e 



aejisc©} 



(2.33) 



- a€{lSC} 



(2.34) 



(2.35) 



We also notice that the variables [x e ,x e ,z, k\, . . . , k n } in the ISC regions can be identified 
with the 3? n ,@ variables in eqs. (2.12) and (2.13). In fact, the fej are integration variables, 
and can be identify with the fcj's in eqs. (2.12) and (2.13). Furthermore, the x@ variables 
are identical to the x @ in eqs. (2.12) and (2.13), since those equations refer to a singular 
configuration, and (as we have remarked earlier) the mapping of eq. (2.18) is the identity 
in the singular region. It follows that the z variables of eqs. (2.14) and (2.33) are identical. 
Hence, from eqs. (2.15) and (2.16), we obtain 



dz 



(2.36) 



(the 1/z factor in the second equation being due to the Jacobian for the change of variables 

x@ — > or@). 

The choice of the counterterms in eq. (2.19) and of the mapping (2.18) should be such 
that the integrals in eqs. (2.31) and (2.33) are easily performed analytically in d dimensions. 
In this way, the C terms contain explicitly the divergences as poles in e. 

We now write eq. (2.17) as 

(O) = J d$ n £O n (3> n ) [B(*„)+Vb(* n )] 

y d$ n+1 {£O n+1 (& n+1 ) 7e(* n+ i)-£[£O n (* n ) c(* n+1 )] Q } 



+ 



+ I /<**n £<?„(*„) C(*„) 

ae{FSC,s} 



+ Yl J d& nt@ £O n (* n ) C (*„,©) 
a ae{isc©} 



+ 



y d* n , e £O n (* n ) £ ffi ,b(*n,e) + y d*„, e £On(*n) e e ,b(*n,e) ■ (2-37) 



Notice that, in the last line, we have substituted £ — > £ for uniformity of notation. This 
is correct, since, as pointed out earlier, in the phase space of the collinear counterterms we 
have x ^ = x@ . 
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It turns out that it is always possible to write 



0©,b(*n,©) + £ CW(* n , @ ) = G @ (* n , @ ) + 5(1 - z) (*„) , (2.38) 

where (<&«,©) is finite in a! = 4 dimensions. 8 The only remaining poles in e are included 
in the last term of eq. (2.38), and have soft origin. It also turns out that in the quantity 



V(* n ) = V b (*„) + 



(*n) + 5f V (*„) + (* n ) 

aG{FSC,s} 



a> —a. 



(2.39) 



all poles in e cancel. With the notation 



[■ 



1 *n=*n 



(2.40) 



we mean that the argument between the brackets is evaluated for values of the phase-space 
variables equal to 3> ra . Notice that the identification <& n = <& n is possible, since $ ra 
refers to the underlying n-body configuration, that must correspond to the Born term. 
Defining now the following abbreviations 

R = £TZ, C {a) = C^C {a \ G & = £G®, B = £B, V = CV , (2.41) 

equation (2.37) becomes 



(O) 



B(* n ) + V(& n ) 

+ ! d* n+1 {O n+1 (* n+1 ) fl(* n+ i) - Y, [°n(*n) C(* n+ i)] a } 

+ y O n (* n ) G e (* n , e ) + y d*„, e 0„(* n ) G e (* n , e ) , (2.42) 

and it is now suited to be integrated numerically, since all the integrals that appear in it 
are finite and can be evaluated in 4 dimensions. 



2.3 Subtraction formalism using the "plus" distributions 

The subtraction method naturally arises when results of NLO computations are expressed 
in terms of distributions in final-state variables. In order to illustrate this issue, we assume 
now for simplicity that there is just one singular region, and, in 4 dimensions, we describe 
the kinematics of the emitted parton (with momentum k) in the centre-of-mass (CM) frame 
of the incoming partons with the following variables 

e = 2fc°/v^, y = cos6, 0, (2.43) 
8 We point out that (/©, although finite, may contain distributions associated with the soft region 2 — ► 1. 



- 12 - 



where s = (k e + k e ) 2 , 9 is the angle of the emitted parton relative to a reference direction 
(typically another parton), and <p is an azimuthal variable around the same reference di- 
rection. The singular regions (soft and collinear) are associated with £ — ► and y — > 1 
respectively. More generally, in d = 4 — 2e dimensions, we can write 



Ll. „l-e 

K 6 -l-2e /i „.2\-e j/- j„. jo<2-2 



2fc (2vr) d - 1 (4vr) d - 



^-^(l-y^c^y^ , (2.44) 



where 



d-3 

dn d " 2 = (sin 0)~ 2e d0 dfi d_3 , y dn d ~ 3 = ^j^j • (2-45) 

then [£ 2 (1 — y) 7£] is regular for £ — > and y — > 1. The phase-space integral of 1Z is infrared 
divergent, and in d = 4 — 2e dimensions (with e < 0) the singular part of the integration is 
proportional to (see eq. (2.44)) 

f dy (1 - j/)" 1 - jf 1 r 1_2e [£ 2 (1 - y) ^] • (2-47) 
In order to deal with the singularities, one uses the expansions 



If we write 



r 1 "* = ~S(0 + Q) - ^ (^J + + Ofc 2 ) , (2-48) 

(1 - y)- 1 - = ~S(1 -y)+ (j^j + O(e) , (2.49) 
with the usual definition of the + prescription 

I df (l) + /K) ^ M 7 M ' (2 ' 50) 

^(') /(>) j' d ,fflhffl. (2 . 52) 

■i V 1 -y/+ y-i 1-1/ 

Inserting eqs. (2.48) and (2.49) into (2.47), and defining y(£,y) = £ 2 (1 — y)7£, for ease of 
notation, we have 



y^y^y) 



y(e,i) 



-NyG) + -'(¥), 

+ £' % r*(?) + (r^) + 9 «' ! ' )+0(e) - (2 ' 53) 
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The first term on the r.h.s. of eq. (2.53), after integration in y and 4>, gives a contribution 
with the same structure as the virtual term, with which it is combined. Since this term 
arises from the <5(£) factor, it can be easily obtained using the eikonal approximation for soft 
emissions in d dimensions. The second term on the r.h.s. of eq. (2.53) is the contribution 
to 1Z proportional to the delta-function term in eq. (2.49), multiplied by the second and 
third term of eq. (2.48). It gives rise to terms that, in the case of final-state singularities, 
can also be integrated in £ and in </>, and yields terms of the same form of the virtual terms, 
with which it is combined. Also for this term it is not necessary to know 1Z in d dimensions, 
since one can use the collinear approximation in the y — ► 1 limit in order to obtain it. In 
the case of initial-state singularities, the same procedure gives terms of the same form of 
the collinear counterterms, with which they are combined. Finally, a term of the form 

remains, where 

H{(?) + (r^) + [ ?2(1 - )K ]}' (2 ' 55) 

Observe that the l/£ factor in front of eq. (2.55) cancels against the phase-space £ factor 
in eq. (2.54), and that [£ 2 (1 — y)1Z] has no singularities at £ = and y = 1, so that 
distributions in £ and y act on a regular function. 

The procedure outlined in this section is fully general. It can be shown that, defining 
R = £71, one can rewrite eq. (2.42) in the form 

(O) = J d& n O n {& n ) [fl(* n ) + F(* n ) 
+ J d& n+1 O n+ i(* n+ i) -R(* n+ l) 

+ J d* n , e O n (* n ) G e (* n , e ) + J d* n , e O n (* n ) G e (*n, e ) • (2-56) 

By handling the + distributions in 7Z according to the prescriptions (2.50), (2.51) and (2.52), 
one automatically generates the real counterterms, provided the variables £ and y appear 
in the phase-space parametrization. If more than one singular region is present, the real 
cross section is decomposed into a sum of terms, each of them having singularities in no 
more than one singular region. For each term, the phase space is parametrized in such a 
way that the variables £ and y, appropriate to that particular singular region, are present. 

The expression of a cross section in terms of distributions has sometimes the advantage 
that the associated projections are not uniquely fixed. In fact, one has the freedom to chose 
the integration variables other than y and £ at one's convenience. This amounts to choosing 
a different projection. 

2.4 Frixione, Kunszt and Signer subtraction 

In this section, we briefly review the FKS general subtraction formalism, proposed in 
refs. [25, 26], including a few modifications that have been introduced recently (see ref. [14]). 
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In FKS one expresses the cross section for the real-emission contribution as a sum of 
terms, each of them having at most one collinear and one soft singularity associated with 
one parton (called the FKS parton). The singular region associated with the final-state 
parton i becoming soft or collinear to one of the incoming partons are labeled by i, while 
those associated with final-state parton i becoming soft or collinear to a final-state parton 
j are labeled by the pair ij. For each singular region one introduces certain non-negative 
functions 9 S of the (n + l)-body phase space such that 

J> + ^<% = 1. (2.57) 

i ij 

We have two options for the range of the indices in the sums: we can let them range from 
1 to (n + 1) (excluding only the i = j possibility in the second sum), or we can assume that 
the Si and Sij are zero (i.e. they are excluded from the sum) if the corresponding regions 
are not singular. For example, if ij refer to a quark and an antiquark of different flavour 
there cannot be a FSC singularity in this region, and we can set S^ = 0. Also, if i is a 
gluon and j is a quark, we may set to zero Sji, since there is no soft singularity associated 
to j becoming soft. Notice that if i and j are both gluons, both terms Sij and Sji appear 
in the sum, since there is a soft singularity for either parton becoming soft. 
The S function have the following properties 

lim ^ + ^^1 =S im , (2.58) 

(2.59) 

imSji + 5u5j m , (2.60) 

(2.61) 

Thus, in a given soft region, (i.e. if parton m is soft), all Si and Sij with i / m vanish, 
and eq. (2.57) is consistent with eq. (2.58). In a given initial-state collinear region, i.e. 
parton m becomes collinear to an initial state parton, only S rn is non-zero and equal to 
one, consistently with eq. (2.57). In a given final-state collinear region, i.e. if partons i and 
j are collinear, only Sij and Sji can differ from zero, their sum being 1, again consistently 
with eq. (2.57). 
We now write 

i ij 

where 

Hi = S.{R , TZij = SijlZ . (2.63) 

The IZi terms give a divergent contribution (i.e. a contribution which has to be subtracted) 
only in the regions in which parton i is soft and/or collinear to one of the initial-state 




9 The notation of refs. [25, 26] has been slightly changed here in order to simplify the discussion. Functions 
Si and Sij of the present paper play the same role as O' ' and G^' '8(kj T — kf T ) of ref. [26] respectively. 
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partons, and the IZij terms are divergent only in the regions in which parton i is soft and/or 
collinear to final-state parton j. Notice that if we have chosen the option of keeping all the 
S functions different from zero, the IZi and IZij functions corresponding to non-singular 
regions are non-zero but finite. 

Equations (2.58)-(2.61) are the only properties of the S functions used in the analytical 
computations of refs. [25, 26]. Their actual functional forms, away from the infrared limits, 
are only relevant to numerical integrations. 

After the (n + l)-body cross section is decomposed as in eq. (2.62), FKS chooses a 
different parametrization of the (n + l)-body phase space for each term, such that one 
can perform the necessary analytical and numerical integrations in an easy way. The key 
variables in the phase-space parametrization associated with Si are the energy of parton 
i (directly related to soft singularities), and the angle between parton i and one of the 
initial-state partons (directly related to initial-state collinear singularities). For Sij, the 
energy of parton i and the angle between parton i and j (related to a final-state collinear 
singularity) are chosen instead. Therefore, there are only two independent functional forms 
for phase spaces in FKS, one for initial- and one for final-state emissions. 

In the following, we will need a further refinement of the FKS decomposition (eq. (2.57)). 
This is because, in FKS, the e and e collinear regions are both singled out by the Si func- 
tions. In POWHEG we will need sometimes to treat the two collinear regions separately. 
We thus introduce the notation 

Si = S?+S®, (2.64) 

with the properties 

Jim Sf = S im , Jim S® = 0, (2.65) 

&m||&@ &m||&@ 

that refine eq. (2.59). Eqs. (2.62) and (2.63) are modified accordingly. 



2.4.1 The S functions 

In the original formulation of the FKS subtraction, the S functions were defined as sets 
of 6 functions. The different contributions to the real cross section, separated out in 
this way, corresponded to a partition of the phase space into non-overlapping regions. In 
the more recent calculation of single-top hadroproduction in MC@NL0 [14], the S functions 
were instead defined as smooth functions. In view of Monte Carlo implementations, step 
functions should be avoided as much as possible, and therefore we consider here the latter 
approach. We introduce the functions d; L and dij, where i,j = 1, ...,n + 1, with the 
following properties 

di = if and only if Ei = or fcj || k @ or ki || k e , ^ 
dij = if and only if Ei = or Ej = or h L || kj , 

where energies and three-momenta are computed in the centre-of-mass frame of the incom- 
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ing partons. A possible definition of the d's is 

di = (jY E i) (l-cos 2 ^) 6 , (2.67) 

dij = (EiE^j a (l - cos b , (2.68) 

where 6{ is the angle between k{ and fc e , 0jj the angle between fcj and kj, s = (k 9 + fc e ) 2 , 
and a, 6 are positive real numbers. Equations (2.67) and (2.68) can be easily expressed in 
terms of invariants using 



which imply 



fce = ^(1,0,0,1), k e = ^(1,0,0,-1), (2.69) 



Ei = 4= (fce + fce)-*i. ( 2 -70) 
V s 

oosfl^l-^, (2.71) 
Ays 

cos% = l-^. (2.72) 

We now introduce the quantity 
and define 



* = (2 - 74) 
^ = ^7-^^rV ( 2 ^) 



Pdij \Ei + E. 



where /i is a function such that 



lim/i(z) = l, lim /i(z) = , h(z) + h(l - z) = 1 . (2.76) 
For example, one can define 10 

h(z) = (1 ~ Z)C . (2.77) 

y ' z c + {l-z) c y 1 

for some positive c. Notice that the h factor is necessary only if one considers both functions 
Sij and Sji (which is strictly necessary only if both i and j are gluons). 
It is manifest that eqs. (2.74) and (2.75) fulfill eqs. (2.58)-(2.61). 



10 In the original FKS formulation h(z) = 0(1/2 — z), which implies that Sij defined in eq. (2.75) vanishes 
if Ej < Ei. With such a choice, a proof was given in ref. [25] that all the infrared singularities are canceled 
by the FKS subtraction. We have checked that, if a smooth form for h is adopted, this proof goes through 
unaltered by replacing 9(z - 1/2) in eq. (4.84) and (4.85) of ref. [25] with h(l - z). 
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As in section 2.4, in order to separate the © and © collinear regions, we modify the 
previous formulae as follows. We introduce 



df = 



^V^Itcos^) 6 , 



instead of di of eq. (2.67). The definition of V in eq. (2.73) becomes 

1 

dkl 



We then define 



5® 



1 



Vdf 



(2.78) 

(2.79) 
(2.80) 



2.4.2 Contributions to the cross section 

We introduce, for the FKS parton i, the following variables 



p - 2 A 



Di = cos 0i , yij = cos Oij , 



(2.81) 



where di is the angle of parton i with the incoming parton ©, and Oij is the angle of parton 
i with parton j. All variables are computed in the centre-of-mass frame of the incoming 
partons. The phase space for the IZi and IZij contributions is written in d 
dimensions as 



4 - 2e 



ra+l 



i=i 



d d ' l ki 



II (2 7r ) d - 1 2fc° 



,1-e 



(4vr 



d-2 



(2.82) 



where y, fi stands for either j/j, S7j or y^, fijj. The transverse angular variables dQ.f~ 2 are 
relative to the collision axis, while d£lf~ 2 are relative to the direction of parton j. The 
singularities for £ — ► 0, y« — > ±1 or — > 1 are treated along the lines of section 2.3. 
The final expression in the FKS formalism results from the cancellation of the infrared 
singularities which emerge in the intermediate steps of the computation. It thus involves 
only non-divergent terms. 

The contributions to the real-emission cross section, in the notation of eq. (2.56), are 



ij , 



(2.83) 



where 



Hi 



1 

Ti 



1 ( 1 



Tlij — £ 

Si 



2 ve 



1 - Vx 



+ 



i + y< 



[(l-yy) & 2 7^] 



(2.84) 
(2.85) 
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If we need to separate the © and © collinear regions, as discussed at the end of sec. 2.4, we 
have 



i ij 

* f 4{(^(iM [(1TBK?Kf1 }' (287) 

The distributions that appear in eqs. (2.84) and (2.85) are defined as follows 

J!« m (i)r£* m ~ i T e '~ G - (2 - 88) 

1 ° iym f i J = } ^ m-mm*-^* . (2 . 89) 



-l 



l=Fy/,5 7-i 1 =F2/ 



The parameters £ c , ^ and <5 must be chosen in the ranges < £ c < 1 and < <5 /j0 < 2. 
The dependence they introduce in the (n + l)-body contribution is exactly compensated by 
the same dependence in the n-body contribution. In the POWHEG framework it is often 
convenient to use the maximal range of integration, and we will thus also use the notation 

with£ c = l, ( - J =(——) with 5 = 2. (2.90) 



The parametrization of the (n + l)-body phase space, appropriate to the integration of IZ-i 
and IZij, can be chosen as the d = 4 version of equation (2.82), as suggested in the original 
FKS papers. This is however not necessary. Any parametrization of the phase space that 
allows a simple handling of the distributions is acceptable. This freedom is exploited in 
the present work, in order to simplify the formulation of the POWHEG method in the case 
of the IZij contributions, where we make a choice of the azimuthal variables different from 
that of eq. (2.82). 

We now consider the soft- virtual term V in eq. (2.56). We define the set of all the 
n + 2 parton labels for an n-body process 

1= {0,G,l,...,n}. (2.91) 

The virtual contribution Vb of eq. (2.37) is given by 11 



(2.92) 



where 



11 We stress that Vb is the contribution to the cross section due to the interference of the virtual amplitude 
with the Born term. It thus includes the corresponding factor of 2. 
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Notice that, in the second sum on the r.h.s. of eq. (2.92), we sum over i ^ j, and thus, since 
(as we show later) Bij is symmetric, each term appears twice in the sum. The definition of 
the finite part V fln depends upon the definition of the normalization factor J\f, for which we 
have adopted the common choice of eq. (2.93), and from the regularization scheme, that 
we assume to be the standard conventional dimensional regularization (CDR). 12 Finally, 
/j, 2 is the renormalization scale, and Q 2 is an (arbitrary) physical scale that is factored out 
from the virtual amplitude in order to make the normalization J\f dimensionless (thus Vfi n 
depends upon /j, 2 and Q 2 ). 

The symbol fi denotes the flavour of parton i, i.e. g for a gluon, q for a quark and q 
for an antiquark. We define 

C g = C A , C q = Cg = C F , (2.94) 

llC A -4T F n f 3 

Ig = g . lq = Iq = ' ( 2 - 95 ) 

, /67 2vr 2 \ „ 23 m , , A3 2vr 2 \ „ 

Vp = ( T - — ) ~ Y T F n f , ,' q = ^ = ( _ - _ j C F . (2.96) 

In case i is a colorless particle all the above quantities are zero. 

The quantities Bij, commonly referred to as the colour-correlated Born amplitudes, 
are defined in the following way 

e « = -h NvmD Xs.s e V m ™ K>)«i T "-< ^ • (2 - 97) 

spins c i^ c j 
colours 

Here -M{ Ck } is the Born amplitude, and {c^} stands for the colour indices of all partons 
in X. The suffix on the parentheses that enclose -M{ Cfc } indicates that the colour indices 

of partons i,j G X are substituted with primed indices in M^y. Furthermore N sym is 
the symmetry factor for identical particles in the final state, D @ are the dimension of the 
colour representations of the incoming partons (3 for quarks and 8 for gluons) , and 5® are 
the number of spin states. The factor l/(2s) is the flux factor. We assume summation 
over repeated colour indexes (c& for k £ I, d i , c'j and a) and spin indices. For gluons T® b = 
if cab, where f a b c are the structure constants of the SU(3) algebra. For incoming quarks 
T^g = t^p, where t are the colour matrices in the fundamental representation (normalized 
as Tr[tt] = 1/2). For antiquarks T£* = — tp a . It follows from colour conservation that Bij 
satisfy 

£ Bij = C fj B. (2.98) 
The soft-virtual term in eq. (2.56) is given by 

V = CV, V = g( QB+J2 ZijBij + VaA . (2.99) 

i+3 



12 If we instead use the dimensional reduction (DR) scheme, we have Va n = Vni R — as/(27r) B 7(/i), 
where, j(g) = N c /6 and y(q) = (iV c 2 — l)/(4iV c ), where iV c = 3 is the number of colours. 
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The quantities Q and Xjj depend on the flavours and momenta of the incoming and outgoing 
partons. They are defined as follows 



l'h ~ lo S 7^2 (ifi ~ 2C fi lo S 



2Ej 



2E. 



+ 2C h lo § -7= ~ lQ g ~ 2 Vi log -± 



2Ei 



~ lo § %p [Ve + 2 C/® l°g& + 7/ e + 2C/ e log£ c 



~" l J 



1 log 2 ^£ + W ^£ log kj ' ki 

2 8 Q 2 ^ 8 Q 2 8 2^ 



Lio 



+ lo S 



2 &j ' &i 



- log 1 - 



kj ' k^ 
2EjE{ 



kj ' k% 
2EjE{ 



kj ■ ki 

log 



2EjEi 



(2.100) 



(2.101) 



where E^ is the energy of parton i in the partonic centre-of-mass frame. 

We finally report the expressions for the initial-state collinear remnants that appear 
in eq. (2.56). For each collinear singular configuration, relevant to the process, we have a 
term G 9 = CQ m (and a corresponding one for G e = CQ e ), where 



2vr 



(1- z)P f ® f ®(z,0) 
~dP f ® f '®(z,e) 



de 



sSj 

i- 24 log « +: 

-K f ® f ®(z) }B f ® fe (z), 



log(l - z) 
1 - z 



(2.102) 



e=0 



for a process in which an incoming parton © of flavour / e splits into a parton f' m (with 
fraction z of the incoming momentum) that enters the Born process B. The superscripts 
on B denote the flavours of the incoming parton, and the z dependence is to remind that 
the incoming © momentum is rescaled. The distributions are defined as in eq. (2.88), with 
£ = 1 — z. The functions P(z, e) are the leading order Altarelli-Parisi splitting functions in 
d = 4 — 2e dimensions, given by 



P™(z, e) = C F 

P q9 (z,e) = C F 

P 9q {z,e) =T F 1 

P"(z,e) = 2C A 



1 + z n y 
l + (l-z) 2 



— ez 



2z(l - z) 

1 - e 
z 1 — z 



l-z 



+ 



z 



+ z(l-z) 



(2.103) 
(2.104) 
(2.105) 
(2.106) 



The distributions K^' control the change of scheme in the evolution of parton distribution 
functions. They are defined in ref. [25], and equivalently, with the notation g in ref. [27]. 
They are identically zero in MS. 
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2.5 Catani and Seymour subtraction 



In this section we briefly review the general subtraction formalism proposed in ref. [27], 
called dipole subtraction. A dipole is defined by three partons: the emitted, the emitter and 
the spectator parton (the last two forming the dipole). In the dipole formulation, one given 
singular region receives, in general, contributions by several dipoles, differing among each 
other by the spectator parton. Thus, the counterterms are associated with dipoles, 
rather than singular regions. The maps M.^ of eq. (2.18) in the dipole formulation (that 
are summarized in section 6) are constructed in such a way that, in most cases, they affect 
only the momenta of the dipole partons, and all other momenta remain unchanged. 13 The 
maps M*") appropriate to the dipole subtraction will be discussed in section 6. These 
definitions, together with the definitions of the relative dipole counterterms (to be found 
in the original reference [27]) are necessary to define a POWHEG implementation. The 
last two missing ingredients are the soft-virtual term V, and the collinear remnants Q @ . In 
this section we report explicitly the form of these terms, expressed in our notation. 
For the soft-virtual contribution V we obtain the following 



V 



2^ 



v «. - E 



— log 

2 * kj ^ ft ^k'l * kj 



B, 



E 



IT 



(2.107) 

Equation (2.107) has been obtained by a suitable manipulation of eqs. (C.27) and (C.28) 
of ref. [27]. The virtual term V fln coincides with that of our eq. (2.92). 
The collinear remnant is given by 



i=i 



l-z 



+ 5{l-z) 



p/e/^l^tfj^log 



i4 



2 X Aim * kf: 



(2.108) 



and the analogous one for Q e . This formula is the translation in our notation of eq. (10.30) 
of ref. [27]. The definition of the functions K, K, Kp.s. and P are given in appendix C of 
ref. [27]. 

In eq. (2.108), £?® e ^ e is the collinear remnant contribution for flavours /©,/ e of the 
incoming partons. Analogously, the superscripts in B (and Bij) single out a given flavour 
combination for the incoming partons in the Born amplitudes and its colour-correlated 
components. 



3 The only exception is when the emitter and the spectator are the two incoming partons. 
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3. NLO with Parton Showers 



3.1 Parton Shower Monte Carlo programs 

A detailed discussion of the ideas upon which Shower Monte Carlo programs are based is 
beyond the scope of the present paper. The interested reader can find some pedagogical 
introductions e.g. in refs. [29, 30]. Here, we simply need to recall few basic features of SMC 
programs. 

First, an MC starts from a kinematic configuration ("hard") which is generated ac- 
cording to an exact LO computation. Usually such configuration is that of a 2 — ► 2 partonic 
process. The final-state multiplicity is then iteratively increased, by letting each initial- and 
final-state parton branch into a couple of partons with a probability related to a Sudakov 
form factor. Thus, if at a given stage of the shower, the scattering process is described by 
m partons, the algorithm decides with a certain probability whether branching is over at 
this stage, or further branchings will take place. In the latter case, one of the m partons 
splits into a pair, generating an m + 1 body final state. Thus, the algorithm defines a 
mapping 



that is fully analogous to the mapping (2.28). Also in this case there is one mapping for 
each singular region, where the singular region is associated with the parton that undergoes 
the splitting. Observe that mappings defined in eq. (3.1) act non-trivially also on the 
momenta of the partons that do not undergo any splitting process. This is due to the fact 
that momentum conservation must be restored after branching, an operation that is usually 
referred to as "momentum reshuffling" in the SMC jargon. Also the value of the momentum 
fraction of the incoming partons may require readjustment, which leads to the fact that the 
value of the luminosity used for the cross-section computation does not correspond exactly 
to what one would have used if the (m + l)-particle matrix element had been computed 
with standard methods. These readjustments usually amount to corrections beyond the 
(leading log) Monte Carlo accuracy, but should be analyzed carefully if one aims at NLO 
accuracy. 

3.2 Including NLO corrections into Monte Carlo programs 

The embedding of a NLO computation into an MC framework, as first clarified in ref. [2], 
aims at reaching NLO accuracy for inclusive observable, maintaining the leading loga- 
rithmic accuracy of the shower approach. This requires that the hardest emission that 
is generated has the correct distribution also far from the collinear directions, and that 
integrated quantities around the soft and collinear directions have NLO accuracy. This re- 
quirements are met in the MC@NL0 approach by carefully tracing the differences of the MC 
simulation relative to the exact NLO one. The similarity of the mappings (3.1) and (2.28) 
are the starting point for this task. The shower algorithm is analyzed to determine its own 
approximate NLO structure in the subtraction framework, in order to determine unam- 
biguously the difference with the exact NLO formulae. 




(3.1) 
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In the POWHEG approach, one performs the generation of the hardest event with 
NLO accuracy, in a framework that does not depend upon the SMC's shower algorithm. 
This is why it is fully independent from the SMC. Furthermore, the subsequent showers 
takes place at softer transverse momenta, and thus affects infrared-safe observables only 
at the next-to-next-to-leading order (NNLO). Thus, the matching problem considerably 
simplifies, since it no longer requires a detailed examination of the properties of the SMC. 

3.3 POWHEG 

In the POWHEG formalism, the generation of the hardest emission is performed first, 
using full NLO accuracy, and using the SMC to generate subsequent radiation. We give 
here a simple illustration of the method, ignoring, for the moment, the complications due 
to the presence of several singular regions in the NLO cross section. We begin by defining 



B(* n ) = B(* n ) + K(*„) 



+ 



d^ d [J2(* n+ i) - C(*„+i)] + J — [G e (* n , e ) + G e (*„, e )] 



(3.2) 



where we have assumed that all the ^n,@ are expressed in terms of the barred 

variables. Next we introduce the Sudakov form factor 14 



A (#„,*) - exp {- / [*^*~0»(M*»->-*>r 



'n — ^ n 



} • (3-3) 



The function k T (<fr n+ i) should be equal, near the singular limit, to the transverse momen- 
tum of the emitted parton relative to the emitting one. The POWHEG cross section for 
the generation of the hardest event is then 

^ = %)# n A($ n ,pf)+A($ n ,fc T ($ n+1 )) ^^W adl , (3.4) 

where it is assumed that <& n +i is parametrized in terms of ^rad and <fr n , and that values 
of k T (<& n+ i) < p™ in are not allowed. The cross section (3.4) has the following properties: 

• At large k T it coincides with the NLO cross section up to NNLO terms. 

• It reproduces correctly the value of infrared safe observables at the NLO. Thus, also 
its integral around the small k T region has NLO accuracy. 

• At small k T it behaves no worse than standard Shower Monte Carlo generators. 

Thus, it fulfills the requirement of the previous subsection for the inclusion of NLO correc- 
tions in an SMC. 



14 Torbj6rn Sjostrand has pointed out to us that a similar Sudakov form factor is also used in PYTHIA for 
weak vector-bosons decay and production, in order to implement a matrix-element matching for the first 



emission in the shower, see refs. [30, 31]. 
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As it stands, the POWHEG formula (3.4) can be used to feed a SMC program, that 
will perform all subsequent (softer) showers and hadronization. If the SMC is ordered in 
p T , we simply require that the shower is started with an upper limit on the scale equal 
to the k T of the POWHEG event. In case the SMC uses a different ordering variable, a 
problem arises, since the POWHEG cross section requires the emissions with higher k T to 
be suppressed in the SMC. This problem typically arises when interfacing POWHEG to 
angular ordered SMC's. It is dealt with by vetoing emissions with larger k T in the shower, 
and by introducing vetoed truncated showers (see ref. [4]), that compensate for the fact 
that in angular ordered shower the hardest emission may not be the first. 

Modern SMC programs, such as HERWIG and PYTHIA, have the capability of generating 
a vetoed shower. This is not the case for the vetoed truncated showers. We point out, 
however, that the need of vetoed truncated showers is not specific to the POWHEG method. 
As discussed in ref. [4] , it also emerges naturally when interfacing standard matrix element 
calculations with parton shower, as in the approach of ref. [32]. At present, there is no 
evidence that the effect of vetoed truncated showers may have any practical important. 



In order to implement the POWHEG method one must specify the separation of the singu- 
lar regions, and the kinematics that associates a given (n + l)-body singular region with an 
n-body one. We discuss POWHEG in the framework of a generic subtraction formalism. 

When one aims at the construction of an event generator, flavour should be carefully 
tracked, since different flavour structures always give rise to different events. We thus 
distinguish the contributions to the cross section also by their flavour structures, which are 
determined by the flavours of the incoming and outgoing partons. We call equivalent two 
flavour structures that differ only by a permutation of final-state partons. In particular, 
we label with the index the flavour structure of the n-body processes and write B^ b and 
V fb for the various Born and soft- virtual contributions. 

We label with a T a particular contribution to the real cross section that is singular in 
only one singular region of integration and has a specific flavour structure. Thus, to each 
a T corresponds one and only one singular region. We then write 



A similar separation also holds for the counterterms, so that they are also labelled by an 
index a r . 

In the FKS case, for example, the a r contributions are obtained by first separating the 
real contribution R into the sum of all its flavour components. For each flavour component, 
one constructs the S functions, according to the procedure of section 2.4.1, and then 
multiplies it by the factors Si or Sij. In the CS case, for each flavour component of the 
real contribution, one defines 



4. The POWHEG method 




(4.1) 



S a — 



(4.2) 
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where a ranges in the set of dipoles T> with the same flavour structure. 

To each contribution a T we can associate an underlying n-body process, with a specific 
flavour structure. The association is performed as follows. If the singular region is collinear, 
the two collinear particles are merged into a single particle in such a way that flavour is 
conserved. In particular, a gg pair is merged into a g, a qq pair is merged into a g, and a 
19 (Q9) P arr is merged into a q (q). If the singular region is soft, the soft gluon is removed. 
Observe that, for non singular limits (for example, in case two quarks, or a quark and an 
antiquark of different flavour become collinear, or in case a quark becomes soft) the flavour 
structure of the underlying n-body process is undefined. 

The factorization remnants also have a flavour structure. We label it with the index 
a@ , and also for these configurations there is an underlying n-body flavour structure. We 
call {a r |/fe} and the set of all values of the indices a r and q@ that have the flavour 

structure of the underlying n-body-configuration equal to fi,. 

We rewrite eq. (2.42) according to the notation of this section (making use of a straight- 
forward extension of the context square brackets introduced in eq. (2.21)) 



<°> = E 



a r e{"r|/b} "eG{ae|/t,} Oi e £{a e \f b } 



(O)g = J d* B O n (* B ) [fi(* n )+V(* n )] , 
(0)g%= /d*n, e O„(* B ) 



n,®) i 



(O) 



R 



n+1 



O n+1 (* n+ i) 22(* n+ i) - O n (*n) C*(*n+l) 
According to ref. [4], we now perform the following manipulation 

(0) a n = (0)% n + (0)l n+1 , 



<o)% n 



(on 



R,n+1 



d$ n+1 O n (* n ) [R (* n+1 ) - C (* n+1 ) ] 

d* n+ ii?(* n+ i) |o„+i(* n+ i) - O n (* r 



, (4.3) 

(4.4) 
(4.5) 
(4.6) 

(4.7) 
(4.8) 

(4.9) 



The term of eq. (4.9) involves real radiation. All other terms (i.e. (4.4), (4.5) and (4.8)), 
have n-body kinematics and, according to ref. [4], should be all lumped together into an 
n-body kinematics term, that was called B. However, we should now carefully distinguish 
the contributions to B according to their flavour structure. We first rewrite eqs. (4.5), (4.8) 
and (4.9) in the following form 



(0) a G %= |d# B B (* B )^G^(* B>e ), 



(0)% n 



J d<& n O n {$ n ) d$ rad {i?(* n+ i)-C(*„ +1 )} 
J d$ n d<S> Tad R($ n+1 ) {0„+i (*„+i) - O n ($ n )} 



(4.10) 
(4.11) 

(4.12) 
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According to ref. [4], we can write the B functions, one for each flavour configuration, as 



B^(*„) = [B (* n ) + V (* n )] A + E / [ d *™d { R " C 

are{a r |/ b } 

+ E + E /f<£ e (*»,e), 



«©e{«<Bl/tl 



«ee{ael/J ' 



so that 



| d* n o n (* n ) bA(*„) = <o>§+ e <o>%„+ 



and 



+ E / d*nd*radi2(*n+l){On+l(*n+l)-On(*n)} 



(4.13) 



^ -/ Ge - W e , (4-14) 

ore{a r |/i,} a®e{o!el/b} "eG{ael/6} 



(4.15) 



We now define the Sudakov form factors 



A /b ($ r 



exp 



arG{a r |/6> 



d$ rad i? (*„+l) 9 (/c T (* n+ i) - pr) 



Bf» (*„) 



(4.16) 

Notice that the identification = $ n is a sensible one only if the underlying n-body- 
process flavour structure of a r is equal to fa. In eq. (4.16), k" T is a function of the kinematics 
variables that depends upon the particular singular region we are considering (its a T index 
is omitted in eq. (4.16) thanks to the context convention). For initial-state collinear sin- 
gularities, we require that k T is proportional to the transverse momentum of the emitted 
parton with respect to the beam axis in the collinear limit. For final-state collinear sin- 
gularities, assuming that the singular region corresponds to momenta fcj and kj becoming 
collinear, we take as k T the (spatial) component of fcj (or equivalently kj) orthogonal to 
the sum ki + kj. In the following we assume that the transverse momentum is computed 
in the CM frame of the colliding partons. 

The factorization and renormalization scales adopted in the definition of B, eq. (4.13), 
and in the definition of the Sudakov form factors, eq. (4.16), are different. In the definition 
of B one adopts a choice that is appropriate to the Born cross section. In the Sudakov 
exponents one must instead adopt a scale of the order of k T . In section 4.4 we show that, 
with the above choice of scales, the Sudakov form factor of eq. (4.16) is equal, at least to 
the leading logarithmic (LL) level, to the DDT [33] Sudakov form factor, and that, in some 
cases, with a simple prescription, one can reach NLL accuracy. 
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The formula for the full POWHEG cross section is 

da = ^BA(*n)#n{A A (* (l)P r) 
fb I 



[d$ rad 6 (k T - p? lin ) AA(* n , k T ) R (* n+1 ) 
+ S ^ « J' < 4 " 17 ) 

a r e{a r |/ 6 } V ™ ; J 

where, for ease of notation, we have dropped the argument in /c" r . The p™ m value 

introduced here is a lower cut-off on the transverse momentum, that is needed in order to 
avoid to reach unphysical values of the strong coupling constant and of the parton-density 
functions. 

As discussed at the beginning of section 2.2, in case the n-body cross section possesses 
singular regions, the observable O n+ \ should vanish fast enough if & n +i approaches two 
singular regions at the same time. Notice that, in the POWHEG cross section given in 
eq. (4.17), the observable function O has disappeared, so that this restriction is no longer 
apparent in the formula. However, thanks to the partition of the different singular regions, 
it is sufficient to apply to the B function a damping factor that suppresses the regions 
where the n-body configuration becomes singular, in order to get a finite result. In this 
way, the POWHEG approach more closely resembles the standard Monte Carlo generators, 
where the hard leading-order matrix element for jet production is appropriately cut off in 
order to get a finite total cross section. 

4.1 Transverse- momentum ordering 

A word of caution has to be said with regard to the separation of the various singular 
contributions in POWHEG, in cases where the underlying n-body cross section also pos- 
sesses singular regions. In standard NLO calculations, these regions are avoided by simply 
requiring that the physical observables one computes should be finite for the n-body term. 
In POWHEG, this requirement is in general not sufficient to guarantee consistency. We 
should also require that the k T of the generated radiation should not be harder than all 
the &: T 's associated with the underlying Born kinematics. More precisely, we should require 
that radiation with k T larger than the smallest k T of the underlying n-body process should 
be suppressed. The separation of R into contributions from the various singular regions 
achieves to some extent this purpose: R ai is suppressed in any singular region different 
from a T . However, consistency with the treatment of soft singularities requires that the 
suppression should be based upon k T , rather than, for example, virtuality. Thus, in the 
FKS case (for example), a most appropriate choice for the df is 

df = (^) 26 2 6 (1 T cos Oif, (4.18) 

dy = (e^:) 2 b (l-cos%)\ (4.19) 

that correspond to the square of the transverse momentum to the power b, rather than the 
form suggested in eqs (2.78) and (2.68). 



* ar — * 
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4.2 NLO accuracy of the POWHEG formula 

In ref. [4], it was argued that the POWHEG formula yields NLO accuracy for infrared- 
finite observables, and that subsequent showering by an SMC program does not spoil this 
conclusion. The second point is a simple consequence of the fact that no radiation with 
k T larger than that generated by POWHEG is allowed in the subsequent shower, and it 
requires no further discussion. We wish instead to address the first point more rigorously. 
In other words, we want to show in detail that formula (4.17), when used to compute 
an infrared-safe observable, yields the correct NLO accuracy. The reader that finds this 
conclusion already obvious can safely skip this section. 

In order to ease the notation, we will drop the 9 (k T — p™ in ) factor in the POWHEG 
formula, always assuming that this factor is present when there is real radiation. 
If we apply formula (4.17) to an infrared-safe observable O, we have 



<0> = E/ d*n^(*n)|A^(* n ,^ in ) O n (* n ) 

fb 



+ E 



d$ rad AA(* n ,fcr) i?(* n +i)O n+ i(* n+ i) 



« r e{a r |/f,} 



fb J 



AA(* n)J ^)+ £ 



[Jd<S> lad AA(* n ,fcr) i*(*„+i; 



a T £{a T \f b } 



Bfb(& r 



O n (*„) 



+ E 

a r e{a r \f b } 



d$ rad A^(* n ,fcr) R(*n+l) (O n+ i(*„ + i) - O n (* n )) 



— ^ n 



(4.20) 



where, in the second equality, we have simply added and subtracted the same term pro- 
portional to R (<fr n+ i) O n (<& n ). We now show that the term in the large squared bracket 
in the third member of eq. (4.20) is equal to 1. In fact 



E 



[7(» rod A-f>(*,„fcr)R(*,, + l) 



a r G{or|/b} 



/ „min 



[y^rad S(k T -p T ) A^(*n,Px) # (*n+l) 



<I> ' — (b 



e{a r |A} 



/ dp T A^(*n,P T )^— J2 



d^rad 0(k T -p T ) R(& n+1 ] 



'Pt 

coo 



O r e{Or|/6} 



/ dpT A f»(* n ,p T ) = 1 - Af»(* n ,pf») , 



(4.21) 
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where we have used the fact that A* 6 (<!>„, oo) = 1. Furthermore, in the last term in the 
large curly bracket of eq. (4.20), small k T values in the integral are suppressed by the 
O n+ i(3> n+ i) — O n (<fr n ) factor, and therefore we can replace A^> — ► 1 and B — > B up to 
higher orders in a s . Equation (4.20) thus reduces to 



<°) = E/ ^n|^ /b (*n) On(*n) 

^$ rad i?(* n+ i) (0„+l(*n+l) - O n (* n )) 



+ E 

Qfr€{Or|/6} 



(4.22) 



up to NNLO corrections. The restriction 6(k T — p" nn ) can now be dropped from the 
<i$ ra d integration, its effect being suppressed by powers of m , and formula eq. (4.22) is 
immediately found to agree with eq. (4.15), thus concluding our proof. 

4.3 Practical implementation of the POWHEG formulae 

The POWHEG cross section (eq. (4.17)) looks very complex, but, in fact, from a numerical 
point of view, it is quite easy to implement using few well-known Monte Carlo techniques. 

4.3.1 Generation of the Born variables 

To begin with, we must generate a Born-like configuration (a point in the <& n space) and 
a value of the index with a probability given by B^ b (^ n ) d<& n . The standard Monte 
Carlo technique used in this cases is the hit-and-miss procedure: one finds an upper bound 
to the cross section, generates randomly the phase-space point, and then accepts it with a 
probability equal to the ratio of the value of the cross section at the given point over the 
upper bound value. This technique is inadequate for our case, since each evaluation of the 
B function requires an integration over the radiation variables. We thus proceed as follows. 
For each singular region, we parametrize the radiation variables $ ra( j in terms of a set of 
three variables in the unit cube, that we call X ra( j = j-X"^, X^ d , -X^j- Similarly, the z 
variable in the collinear remnants is parametrized in terms of one of these three variables, 
that we take to be X^ d . We then introduce the function 



B A (*„,X rad ) = [B (<&„) + V (*n)] f b 

d$rad 



+ E 

ar€{a T \f b } 



+ 



E \ 



5-Xrad 

dz 



{R(* n+1 )-C($ n+1 )} 



a©€{a©|/i,} 



(i) 

rad 



(* n>e ) + E 7 



a e G{« e l/i)} 



dz 



dX, 



(i) 

rad 



G% e (*„, e ), 
(4.23) 



so that 



SA(* B ) = f 1 dX ^ d f dX™ f dX^ d BA(* B> X rad ) , (4.24) 
Jo Jo Jo 
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and we define 



(4.25) 



fb 



There are computer programs that, after performing a single integration on a given func- 
tion, can efficiently generate points in the integration range, distributed according to the in- 
tegrand function. One such popular program is the BASES/SPRING package [34]. The adap- 
tive Monte Carlo integration routine BASES performs the integration of the non-negative 
function, and stores the necessary intermediate results. The routine SPRING uses these 
information to generate unweighted events. In our case, one integrates the B(<f> n , -X" r ad) 
function in the full (<& n > -^rad) space using BASES. Then one generates (3? n , -X'rad) points 
using SPRING. For each generated phase-space point, one chooses an fb value with a prob- 
ability equal to Bf b (<& n , X xa a)/B(<& n , X ra( j). At this point, the X ra( j values are discarded, 
and one has generated the (<& n , fb) values with probability proportional to B^ b {<^ n ). In 
essence, by doing a single (n+ l)-body phase-space integration, one is able to generate the 
Born configuration with reasonable efficiency. 

As already pointed out, if the Born configuration possesses singular regions, the n- 
body phase space should be suitably constrained, in order to avoid them. This is also 
the case in standard SMC programs: when one deals with cross sections with singular 
matrix elements, typically cross sections that include the production of a light parton, 
one specifies a transverse-momentum cut on its momentum, in order to get a finite total 
cross section. Alternatively, a weight W($* n ) should be attached to the B function that 
suppresses the n-body singular regions, so that the integral of W x B is finite. Events are 
then generated using W x B instead of B, and a weight VF _1 (<& n ) should be attached to 
each event. For example, in order to generate a sample of Z + jet events, knowing that we 
will select jets with transverse energy greater than EjP, we can restrict the Born events so 
that Et > -E^ ut /2, Et being the transverse energy of the radiated parton. Alternatively, 
we weight B with E T , attach the weight 1/E T to the generated event, and unweight the 
events after the event sample (including cuts) is generated. The cut on the transverse 
energy of the associated jet should effectively cut off the events with low parton Et, so 
that the unweighting will be really possible. 

4.3.2 Generation of the hardest-radiation variables 

Given the Born kinematics (<fr n , fb), we must now generate the hardest-radiation configu- 
ration, characterized by (a r , 3>^d)' w ith a r G {a r |/b}, with probability 




71 



,£*r 

rad ' 



(4.26) 



The Sudakov form factor can be written as 



A A (*n,Pr)= ft A «r(*n,P T ), 



(4.27) 



a I e{a I \f b } 
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where 



A&(*n,Pr) =expi - 



-Pt) 



(4.28) 



Under these conditions, the problem of generating the radiation variables according to 
eq. (4.26) can be reduced to the problem of generating them with probabilities 



rad ' 



(4.29) 



by using the highest-bid method, illustrated in appendix B. We are thus left with the 
problem of generating radiation variables according to eq. (4.29) for a fixed value of a T . 
This problem can be dealt with using the veto technique, illustrated in appendix A. In 
order to use this technique, we need a sufficiently simple upper bounding function 



n+lj 



<F(** * B ). 



(4.30) 



BA(* n ) 

This can be found by taking the singular limit of the left hand side of eq. (4.30), that has, 
in general, a form suggested by the factorization theorem, and by elementary properties of 
the parton densities in the case of initial-state singular regions. Once the functional form 
of F is guessed, its normalization is found by scanning the 3>„+i phase space. 

Within a given subtraction method, one has typically two kinds of upper bounding 
functions: one for final-state radiation and one for initial-state radiation. In section 7 
we illustrate explicit forms for F in the FKS and CS frameworks, for both initial- and 
final-state radiation. 

Notice that, in general, it is not necessary to separate out all a r regions in order to apply 
the veto method. In many cases several regions can be group together, thus simplifying 
the generation algorithm (see section 7). 

4.4 Sudakov form factors and NLL soft gluon resummation 

The purpose of the POWHEG method is to reach NLO accuracy for inclusive quantities, 
and leading logarithmic (LL) accuracy for exclusive final states. In this section we ad- 
dress the following question: to what extent we can do better than LL accuracy in the 
POWHEG framework? First of all, the POWHEG method deals with the hardest emission 
only. Subsequent emissions are handled by the shower Monte Carlo to which POWHEG is 
interfaced, and will have, in general, only LL accuracy. However, exclusive observables that 
are especially sensitive to the hardest emission will benefit from an improved POWHEG 
accuracy. 

In this section, we show how to improve the logarithmic accuracy of POWHEG. 
We show that in many cases this improvement requires very minor adjustments of the 
POWHEG Sudakov form factor, and that in general, NLL accuracy at least in the large 
N c limit (where N c is the number of colours) should be easy to achieve. In all this section, 
for ease of notation, k will denote the momentum of the soft gluon. 

The results of this section can be summarized as follows: 
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1. The POWHEG Sudakov form factor is accurate at the LL level, provided that the 
strong coupling constant and the parton density functions in the Sudakov exponent 
are evaluated at a scale of order k%. 

2. In case of processes involving no more than 3 coloured partons (in the initial or final 
state), NLL accuracy is achieved by replacing the strong coupling constant in the 
Sudakov exponent with 



where the MS, 1-loop expression of a s should be used. Furthermore, the parton 
densities in the exponent must be evaluated at a scale of order k%. The argument of 
a s in eq. (4.31) can also be taken equal to a function of the radiation variable that is 
of order k\ in the soft or collinear limit, but becomes exactly equal to k% in the soft 
and collinear region. 

3. In case of processes involving more than 3 coloured partons, the procedure of item 2 is 
not sufficient to guarantee NLL accuracy. There are in fact soft (NLL) contributions 
that do exponentiate only in a matrix sense, so that, in order to deal with them using 
standard Monte Carlo techniques suited for the evaluation of ordinary exponential 
(like, for example, the veto method), one should diagonalize their colour structure. 

4. We will show that, for processes involving more than 3 coloured partons, the correct 
exponentiation of the soft (NLL) contributions discussed in item 3 can be easily 
recovered for the dominant terms in the large- N c limit. 

In the rest of this section we will demonstrate the above points. We assume, in the following, 
that the reader is familiar with Sudakov resummation techniques, and in particular with 
ref. [35]. We also reassure the reader that, if she/he is willing to accept the above points, 
and is not interested in implementing the 4th point of the above list, the rest of this section 
can be safely skipped. 

We begin by recalling the structure of Sudakov resummation of soft gluon effects in 
QCD, taking the notation and the results given in ref. [35]. One usually assumes that there 
are kinematic constraints that limit the emission of large transverse-momentum partons, 
called Sudakov weights u in [35]. In the POWHEG Sudakov form factor, the constraint is 
given by the requirement that all radiation processes have transverse momentum less than 
a given p T , so that u = #(p^ ~~ &t) m our case. We now review the ingredients that build 
up the Sudakov resummation factors at NLL level. Radiation from a final-state massless 
parton i carries a factor Ji(p T ), given in eq. (10) of ref. [35] 





(4.31) 




where the function A is defined in eq. (4.31). Furthermore, 



k ■ n 



(4.33) 



z = l- 



ki ■ n 
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where n is a timelike vector, that we take to coincide with the time direction in the partonic 
CM system, 15 and k T is the transverse momentum of k relative to hi, defined as 

k\ = 2(1 - z)h ■ k (4.34) 

in [35]. It corresponds in the soft-collinear limit to the transverse momentum of k with 
respect to h L in the CM frame of the incoming partons. Pj are appropriate combination of 
the Altarelli-Parisi splitting functions, defined in eq. (11) of ref. [35]. The soft (i.e. z — ► 1) 
and collinear singularities in formula (4.32) give rise to logarithmic terms, with a structure 
(at the NLL level) 

oo oo 

log MP?) = E c f L «s(M 2 ) log fe+1 ^ + E^W) log* ^, (4.35) 
k=i V k=i V 

where \i is a reference scale of the order of the renormalization scale. 

The logarithms arise in the following way. The collinear region of integration generates 
a logp T . The soft region (i.e. z — > 1) also generates a logp T , which arises from the 1/(1 — z) 
terms in the Pi(z) functions. The NLL expansion of a s (fc 2 ) in powers of a s evaluated at a 
reference scale \x has the structure 



a s (kl) = a s (/J 2 ) 



E c ? LL («s (^ 2 ) log ^) + E c " ,NLL ° s (^ 2 ) ( as (v 2 ) log ^) + . . . 



(4.36) 

and thus generates higher powers of logarithms. The coefficients c^' LL in eq. (4.35) depend 
upon the coefficient of the 1/(1 — z) term in Pj, and upon the c"' LL terms in eq. (4.36). The 
coefficients c^' NLL depend also upon the full form of Pj, since non-soft, collinear terms can 
generate single-log contributions, and upon the C J' NLL terms in eq. (4.36). Furthermore, 
the 0(a 2 ) term in A also contribute to the c J '^ LL coefficients. It arises from the z — > 1 
singular part of the Altarelli Parisi splitting function, that has the form [36] 



67 ^ 2 . 

C A nf 

18 6/9 



c g = 2C A , c q = C F . (4.37) 



This relation holds for both the spacelike and timelike splitting function. The correspond- 
ing a 2 correction is associated with a collinear and soft singularity, and thus yields two 
logarithms. When combined with the c"' terms of eq. (4.36) it gives rise to NLL con- 
tributions of order a^log 2 p T /fi, ccf log 3 p T /fi, and so on. We notice that the expression in 
eq. (4.32) can be obtained from its own C(a s ( / u 2 )) expansion 

logJ^ T ) = -4vr J ^S + (k 2 ) e(kl-p 2 T ) ^«s(m 2 ) Pi(z) + O (ai( M 2 )) (4.38) 

provided one replaces a s (fj> 2 ) — > A(a s (k 2 )). 

We now notice that the POWHEG Sudakov form factor has an exact expression of order 
a s in the exponent. Since formula (4.38) generates the dominant log 2 p T and \ogp T terms, 



3 In ref. [35] one n is defined for each ki, Here, for simplicity, we identify all of them. 
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it must also be contained already in the POWHEG Sudakov exponent. It thus follows that, 
with the replacement of eq. (4.31) in the POWHEG Sudakov, we automatically include all 
LL and NLL terms of the J(p T ) factors. A similar reasoning holds for the resummation of 
initial state emissions 16 , called A factors in [35], and leads to the conclusion that, besides 
the replacement (4.31) one also needs to evaluate the parton densities in the Sudakov 
exponent at a scale of order k%, in order to achieve NLL accuracy. 

The exponentiation of genuine soft, non-collinear interference contributions is more 
problematic. These contributions do not generally factorize in terms of the Born cross 
section. The real emission cross section in the soft limit has the well known structure 



1Z « 4-7TQ s 



5>. 



kj ' k -j 



(y kj * k^j [kj * k ^ 



+*E 



(h ■ k)' 



Ci 



(4.39) 



where Bij are the colour correlated Born amplitude, defined in eq. (2.97), and Cj is the 
Casimir invariant of the colour representation of parton i. If parton i is massless, also 
collinear singularities are present in formula (4.39). They can be separated out by exploiting 
the techniques used to derive the results in ref. [35]. Since that technique is illustrated in 
unpublished notes, in the following we sketch its main points. Defining k{ + kj 
have 



k 



we 



kj ' k -j 



ki ' k 



+ 



kj * kij 



( kj ' k ) ( kj m k ) [ kj ' k ) [ kj j * k ^ [ kj * k ^ ( kj j * k ^) 

In case of massless partons, for example, kf = 0, we can separate out the collinear compo- 
nent 



kj * k 



1 



ki ' kij 



ki 



n 



k • n 



+ 



h • n 



(ki • k) k • n ' 



(4.41) 



{ki • k)(kij ■ fc) (ki ■ k) _ (kij • A;) 

where n is an arbitrary timelike vector. In the first term of the right hand side of eq. (4.41) 
there are no collinear singularities, since the factor in the square bracket vanishes when k 
becomes collinear to k^. The second term in (4.41) is collinear, but being independent of j 
gives a contribution of the form 



1 kj ■ n 



ocB- 



1 h ■ n 



(ki ■ k) k ■ n (ki ■ k) k ■ n 



(4.42) 



factorized in terms of the Born cross section. Soft terms that factorize in terms of the Born 
amplitude are automatically included in the POWHEG framework, since the POWHEG 
exponent contains precisely the factor 1Z/B. Not so for the interference terms. We are thus 
left with terms of the form 

2(ki • kij kf) 



kj * k> 



ki ■ n 



[kj * k^j ^J^ij ' ^0 k ' tl 



(ki • k) (kij ■ /c) 



(4.43) 



We distinguish in the sum the massless and the massive case. The 6(k T — p T ) function in 
the POWHEG Sudakov exponent, eq. (4.16), could in principle be different for the various 



3 In ref. [22] this case is discussed in details, in the framework of Z pair production in hadronic collisions. 
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collinear regions of the soft parton. However, since formula (4.43) does not carry collinear 
singularities, the dominant integration region does not require any small angles. Under 
this condition, k T is of the order of k°, and at the single-logarithmic level one can replace 



e{kl- P %) ^e(k 2 - P 2 T ) 



(4.44) 



In other words, the integral in the radiation variables of the generic term of formula (4.43), 
with a theta function constraint 9(k T —p T ), and k T defined relative to one generic ki yields 
a result of the form alog£> T + b. If we replace the 6(k T — p T ) with #(fc — p T ) the result will 
become alogp T + b' , i.e. the logarithmically enhanced term remains the same. 

The angular integration of the /c-dependent coefficients in eq. (4.43) yields (see Ap- 
pendix C) 



1 



1 



(ki ■ k) 



ki ■ k 



k,; 



n 



2vr (kjj ■ kj) 2 n 2 
kl g kUn-hf 



if k 2 = 0, and 



(4.45) 



k 2 



2ir 



(ki ■ k) (kij • k) 



A 



log 



i + A j 
i-A/ 



(ki ■ kijY 



(4.46) 



if k 2 / 0. The NLL resummed Sudakov form factor associated with these soft emissions 
has the form 



-lint 



exp 



v ^ dk° 



>(k°-p T ) 



(k 2 o) 



4n 



M 



(4.47) 



where S mt corresponds to the expression for E in eq. (9) of [35] (see also eqs. (14) and (15) 
there), excluding the J factors. 17 Here M is the Born matrix element, viewed as a complex 
vector in all its colour components. More precisely, M, is a tensor in colour space, carrying 
a colour index for each parton entering or leaving the amplitude. It thus spans a linear 
space, equipped with a sesquilinear product. The matrices Tf can be seen as operators 
acting on this linear space. In facts, they act on the colour index of the ith particle in 
M.. Thus the exponential is to be seen as the exponential of an operator in colour space. 
In ref. [35] the exponential is written as an energy-ordered one, to remind that this kind 
of exponentiation takes place because leading-log soft emission is dominated by energy- 
ordered graphs, where softer emission take place later (i.e. closer to the final-state lines). 
In fact, at NLL accuracy the ordering has no effect. 

In the POWHEG formalism, the exponentiation of the soft interference terms has the 
form 



S mt = exp 



. 2 / v e (k » _ Pt) M * s. n ^ (a, + w m 

J ■ , ■ k 



47T 



\M\* 



(4.48) 



7 In ref. [35], only the case of massless partons is considered. 
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since in POWHEG the ratio R/B appears in the Sudakov exponent. This does not agree 
in general with eq. (4.47). There is, however, one important case in which it agrees, that 
is to say, when M. is an eigenstate of all the operators Yla^i^j- ^ n this case we can 
replace ^ a T? in the exponentials in eqs. (4.47) and (4.48) with their eigenvalues, and 
the two expressions become identical. It is also apparent that the Bij are in this case all 
proportional to B, and again we have complete factorization of the soft contribution in 
terms of the Born cross section. As discussed in ref. [35], if there are no more than 3 
coloured partons entering or leaving the Born amplitude, it turns out that A4 is always 
an eigenstate of the operator Y^a^T®, f° r an y i>3- We conclude that, in this case, the 
prescription of eq. (4.31) is sufficient to guarantee NLL accuracy. 

If there are more than 3 coloured partons in the amplitude, the standard POWHEG 
formula does not allow one to simply obtain full NLL accuracy. It should be very simple, 
however, to modify it in such a way that the interference soft terms are correctly resummed 
at NLL accuracy, at least as far as the dominant terms in the large- iV c limit are concerned. 
In the large- A^ c limit, the Born amplitude can be written as the sum of independent planar 




Figure 1: The two incquivalcnt planar colour configurations in e + e — > qqgsgA- Notice that the 
square of each amplitude is not symmetric in the exchange of the final state gluons. 

colour structures, that differ only by permutations of the external lines, as illustrated in 
fig. 1. We write 

p p 



where the subscript pi stands for planar, and p labels the different planar colour compo- 
nents. Observe that interference terms in the square of Ai p \ do not appear, since they 
give rise to non-planar structures, and are suppressed in the large N c limit. Soft emission 
factors give rise to leading N c contributions only if they act upon adjacent colour lines in 
the planar structures . In this case their effect amounts to a factor of N c /2, and they do 
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Figure 2: Soft emission colour factor T§T^ does not alter the colour structure of the Born planar 
amplitude, and it acts as a multiplicative constant. The close colour loop provides an extra factor 
oiNr. 



not alter the colour connections of the Born planar amplitudes. It follows that 



cxp 



I-Mpll 2 

i 



P i 2 exp 



M 



pi 



E 

i=j±i 



4vr 2 



>p r pi 



p | 2 exp 



pi 



47T 



Kll 2 



(4.50) 



where with the notation z = j±l we (somewhat imprecisely) indicate that we only consider 
adjacent colour lines for a given planar ordered amplitude. The last line of eq. (4.50) is 
similar to the form that appears in eq. (4.48), except that there is one term for each 
planar amplitude. This suggests how to implement it in the POWHEG framework. In 
section 4 we have introduced a classification of the Born contributions in terms of their 
flavour structure and a classification of the real amplitude contributions in terms of 
their singularity regions and flavour structures a r . We need to extend this classification 
to include also the colour structure in the large N c limit. One simple way to do so is to 
compute B p i and R v \, the Born and real terms in the planar limit, for each of their given 
flavour structure. Then we define 



gfb,p 

fifb,P — £>fb P 1 



fb.p 
P "pl 



JjCt r ,p r _ J^OL r 



Er,a r ,pr i 
pr ""pi 



so that 



Bf» = V B fb ' p , R ar = V R ar > Pr . 



(4.51) 



(4.52) 



Pr 
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The labels p and p r refer to the contribution of a given planar structure in the Born 
and Real contributions respectively. Next we must associate with each (a r ,p r ) pair an 
underlying Born (fb, p) pair. As far as /& is concerned, the association is performed along 
the lines explained in section 4. For p and p r , one proceeds as follows. If the singular 
region is soft, the Born planar colour structure is obtained by simply removing the soft 
gluon (and joining its colour lines) from the real emission colour structure. If the singular 
region is collinear, we need to show first that each R ar 'P r has collinear singularities only 
for collinear particles that are nearby in its planar structure. This property follows from 
the fact that the ratio 

(4.53) 



Pr U p\ 

has no collinear singularities, since the denominator is the large iV c limit of the numerator, 
and therefore it has the same singular structure. It follows then, from the second equality 
in (4.51), that R ar >P<~ has the same collinear singularities of R^' pr . Thus, since planar 
cross sections have collinear singularities only when nearby partons are collinear, R ar >P<~ 
has the same property. Joining the planar colour of nearby partons yields a valid Born 
planar colour configuration. We now introduce 



B 



h,p 



B ] 



£>fb,P 

re! 

ET?fb,p 
p n pl 



(4.54) 



and the corresponding POWHEG Sudakov form factor 



A^(* n)Pr ) =exp 



E 

a r ,p r £{a r ,p r \f b ,p} 



a r ,p r 



Bf>>P(* n ) 



>, (4-55) 



that should correctly generate large angle soft radiation in the large iV c limit. Equa- 
tions (4.54) and (4.55) certainly satisfy the formula (4.50) for the soft interference terms 
in the large N c limit. We must still show that collinear singularities are treated correctly 
for all N c , i.e. that formula (4.55) is equivalent to formula (4.16) in the collinear regions. 
This is the case if „ 

. T>a r ,p r T>a r 

E W = m (4 ' 56) 

p r e{a r ,Pr\fb,p} 



in the a r collinear region. Using eqs. (4.51) we write 



E 



R c 



pr£{a r ,pr\fb,P} 



Bfb'P 





l^p r e{a r ,p r \f b ,p} JX p\ 











(4.57) 



The second term on the r.h.s. of eq. (4.57) yields 1 in the collinear limit. In fact, it is 
easy to convince oneself that, in this limit, both its numerator and denominator yield the 
(planar) Altarelli Parisi splitting kernel, that simplifies in the ratio. Thus A^ f " p ($ n ,p T ) 
becomes p independent in this limit, and the overall B^^p factor can be summed in p, to 
yield back B^. 
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We conclude this section by reminding the reader that the importance of the inclu- 
sion of soft-interference terms should not be overemphasized at this stage. After all, the 
POWHEG formula in eq. (4.48) differs from formula (4.47) only starting at order a^. It 
is thus unlikely that these terms will have important effects in a full POWHEG imple- 
mentation. Nevertheless, if one wants to assess their importance, one can, to begin with, 
implement their large N c resummation and study its impact. 

4.5 Interfacing POWHEG to an SMC 

The POWHEG algorithm generates the kinematics and flavour configuration of the hardest- 
emission event. The event should be fed into a SMC using the Les Houches Interface for 
User Processes [16] (LHIUP from now on). In particular, one should require that no events 
harder than the one generated by POWHEG should be generated by the SMC. This is 
achieved by setting the variable SCALUP of the LHIUP equal to the k T of the POWHEG 
event. The LHIUP specifies how to pass the kinematics and flavour structure of the hard 
event to the SMC. 

4.5.1 Colour assignment 

The LHIUP also requires that the colour connections of the hard event (in the large iV c 
limit) should also be specified. POWHEG does not, in general, generate these large-iV c 
colour structures. They are needed (and can be generated) only if one wishes to reach large- 
N c NLL accuracy of the Sudakov form factor, in events with more than 3 coloured partons at 
the Born level, as discussed in section 4.4. If this is not the case, the generation of the colour 
configuration should be performed after the POWHEG event has been generated. Here 
we illustrate two acceptable ways to perform colour assignment. The first (and simpler) 
approach is the following: 

• Generate the POWHEG event in the standard way. 

• Compute the different (planar) colour contributions to the Born cross section, at the 
kinematics of the generated underlying Born configuration. 

• Pick an underlying Born colour configuration, with a probability proportional to the 
corresponding (planar) colour Born contribution. 

• If no radiation has been generated, this is the colour structure of the event. 

• If radiation has been generated, POWHEG has also generated a a T index, specifying 
a singular region. In this case we always assume that the emitted parton is (planar) 
colour-connected to the emitter. This fully specifies the planar colour structure of 
the event. 

This method only requires the calculation of the planar colour-structures of the Born term. 
In the second method (which is usually implemented in MCONLO) one computes all the 
planar colour contributions to R, and chooses the colour configuration with a probability 
proportional to the corresponding contribution. 
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5. POWHEG in the Frixione, Kunszt and Signer framework 



In this section we construct a mapping from the phase space to the barred and radi- 

ation variables of eq. (2.28) (and its corresponding inverse mapping) for the two kinds of 
singular regions i (initial state) and ij (final state), that is compatible with the FKS sub- 
traction method. These is the only missing ingredients one needs to construct a POWHEG 
generator in the FKS framework. 

We stress that the mapping that we propose is not the only possible one. In practical 
examples one may find convenient to depart from this approach. It is however fully general, 
and as such it may be implemented once and for all in a computer code to be used for all 
processes. 

5.1 Initial-state singularity 

5.1.1 Radiation and barred variables 

We first show how to construct the underlying Born (i.e. the barred variables) and the three 
radiation variables from the 3>„+i kinematics, for the case of the initial-state singularity 
region. 

Without loss of generality, we assume that the FKS parton is the (n + l)th one, so that 
the mapping from ki to ki (i = 1, . . . , n) does not requires a relabeling of the momenta. In 
the FKS formulation, k n+ \ is a function of £, y = cosO and <j) (see eq. (2.43)), so that, in 
the centre of mass of the colliding partons, i.e. the CM of the (x® K® + x e K e ) system, we 
have 

kn+i = ~2~£' k n+ i = fc° +1 (l,sin0 sin0, sin 9 cos<^>, cos#) . (5.1) 

We take £, y and 4> as the radiation variables for the singular region, and obtain 



d 3 k n+1 



2fc0 +1 (2vr)3 (4tt) 

We introduce the momentum 



■ididyd4>. (5.2) 



n 



ktot — ] ki — x®.K® + XqKq fen+i i (5-3) 
i=i 

and construct a longitudinal boost Ml (longitudinal with respect to the incoming beams) 
such that B^fctot has zero longitudinal component. Notice that is unique, being given 
by a longitudinal boost with boost angle equal to minus the rapidity of k to t- Then we 
construct a transverse boost By such that B^Bi/c tot has zero transverse momentum. We 
then define the barred momenta as 

h = B L 1 B T B L h, i = l,...,n, (5.4) 

and define 

n 

hot =Y,ki= ^L^T^Lhot ■ (5.5) 
i=l 
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We notice that, by construction, fctot and fctot have the same invariant mass and rapidity. 
We now define x ffi and x e in such a way that 

x ffi iT ffi + x e K e = k t ot • (5.6) 

An explicit expression for x B and x e is easily obtained, by using the fact that k tot and k to t 
have the same invariant mass and rapidity. We get 



Xffi — x ( 



r— 2-gi + y) _ r— 2-^1-y) 



and 

dx m dx e = ^_ ® . (5.8) 

Observe that we always have x e < x e and x e < x e . The phase space d& n+ \ can be 
rewritten as follows 



(n+l \ n+l 
+ x e K e - ^2 h j n 
i=l / 1=1 



ci 3 ^ 



2*?(2tt) e 



= cfe© cfe e t-^t d£ dy d<p 

(4vr) J 1 - £ 



n \ n 



{2^5 A [x m K S) + x e K e -k n+1 -Y J h \ [J " 1 



, 3 

i=l / i- 



\ 2*?(2T) a 



dx ® dx e "7^3 d ^ d V ^ ( 27r ) 4 <^ 4 ^e^e + x e K e ~Yj ki J II 2 p(2 



7T) 3 



(4vr) 3 1 - £ 



d(,dyd<t>d$ n , (5.9) 



where we have used the boost invariance of the n-body phase space. Thus, from eq. (2.29) 
we obtain 

s £ 

d$ rad = -7^3 ^ ^ ^ ' ( 5 ' 10 ) 

5.1.2 Inverse construction 

We describe now the construction of the full (n + l)-particle phase space, given the barred 
variables # ra and the radiation variables £, y = cos 9 and <f>. Using eq. (5.1) we fully 
construct k n+ ±. Inverting eq. (5.7), we have 



*e / 2-g(l-y) x e / 2-g(l+y) 

9 " p-ai + yV 9 " V 2 - «i - ») ' ( J 

The range for the variable £ is now restricted by the requirement that both x e and x e be 
less than 1. This gives 

max 

(5.12) 
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with 

£max = 1 - max 



2(1 + y) g 

V(i + 4) 2 ( 1 -y) 2 + 1 6y^ e + (1 - y)(l - x%) ' 

2(1 -y) 4 

V(l + x|) 2 (l + y) 2 - 16y i| + (1 + y)(l - x%) 



(5.13) 



From k n+ i and x @ we can construct fc to t = a^e-^e + x e^e ~ k n +i, and summing the n 
barred momenta, according to eq. (5.5), we can compute 

n 

hot = ^2^- ( 5 - 14 ) 

i=l 

The four vectors fctot and ktot have the same invariant mass and rapidity by construction, 
since the relation between the x@ and the x@ was obtained precisely from these conditions. 
We then construct the boost Ml such that !>Lfc tot has zero rapidity. We will also have that 
B^fctot has zero rapidity. Then we compute the transverse boost Mt such that 

B T B L fc to t = M L hot- (5.15) 

Finally, the momenta ki, i = 1, . . . , n are obtained as 

k i = Bl 1 n T 1 B L k i , i = l,...,n. (5.16) 

This completes the construction of the (n+l)-body phase space, starting from an underlying 
Born configuration and the three radiation variables. 

5.2 Final-state singularity 

5.2.1 Radiation and barred variables 

In this section, we show how to construct the underlying Born and the three radiation 
variables, given the 3? n +i variables, in the case of final-state singularity. We assume, 
without loss of generality, that the singular region is associated with partons (n + 1) and 
n, that is, the (n + l)th parton becoming collinear to the nth parton, or becoming soft. 
Our mapping is constructed in such a way that 

= x <b ) = x e • (5-17) 

Thus, the partonic CM frame of the final-state (n + l)-particle system coincides with the 
CM frame of the n barred momenta ki, and we work in this frame from now on. In this 
section (and only here), we need to introduce a short-hand notation for the modulus of the 
space component of a four-vector in the CM frame 

p=\p\. (5.18) 

We then define 

n+l 

q = K m x m + KqXq = , (5.19) 



i=i 
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and, in the CM frame, we have 



q = 0, q 2 = (q ) 2 . (5.20) 



Introducing the four momentum 

k = k n + k n +i , (5-21) 

we define the radiation variables as 

e = ^%i, V = -r^—r^ ' cj> = 4>Uxk,k n+1 xk), (5.22) 

where f/ is an arbitrary direction that serves as origin for the azimuthal angle of k n+ \ 
around k, and "x" is the cross vector product. The notation 4>(vi,V2) denotes the angle 
between v\ and v^- Thus, (ft is the azimuth of the vector k n+ \ around the direction k. 
Notice that only £ and y correspond exactly to FKS variables, since the FKS azimuth is 
usually defined with respect to k n rather than k. We prefer the choice (5.22) because we 
want to defined the mapping in such a way that the k direction (rather than the k n one) 
is preserved. Our choice, however, makes only irrelevant differences in the FKS formalism, 
since the real-emission cross section has singular distributions only in the variables £ and 

y- 

We introduce the recoil four-momentum and mass 

n-l 

krec = ki , M^ cc = ^rec ■ (5.23) 

i=l 

We have 

^rec = q k , k rec = k . (5.24) 

We construct a boost B along the k rec direction, such that the 4-momentum (q — M k rcc ) is 
light-like, that is to say 

(q-Mk rcc ) 2 = 0. (5.25) 
The boost velocity is easily computed 

q 2 + (fef? cc + A; rcc ) 

Since q° = k° + /cj? ec , and k° > k = fe rec , (3 is positive and smaller than one. Thus, the 
boost B always exists. We define the barred momenta as 

ki=Mki, i = l,...,n-l, k n = q - B k rec , (5-27) 

so that they clearly satisfy momentum conservation 



Y,~ k i = 1- ( 5 - 28 ) 
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In order to obtain c?<I> ra( j, we write the (n + l)-body phase space as 



n+i ,31.. / "+ 1 



1 2k° (2vr)3 
d 3 k n+ \ d 3 k 



2k° n+l (2tt)3 2fc0 (2vr)3 ±± 2 *° (2tt) 



rt— 1 

n 



d 3 L 



n-l 



(2vr) 4 5 4 U-fc- J> , (5.29) 



i=l 



where, in the last equality, using eq. (5.21), we have traded k n for k as independent variable. 
We thus have now 



k,, 



k° — k® , i + k® . 



(5.30) 



n+l ' 

We identify the phase space in eq. (5.29) with the phase space written in terms of the 
barred and radiation variables, multiplied by a Jacobian factor J 



d 3 k 



n+l 



n+l 



d 3 k 



n—1 



2k» +1 (2„)3 2k« (2nf y 2k? {2*f 



n 5i ££s(W«-*-i> 



= d$ rad eZ$ n = (Jd^dydcp) 



n 



i=i 



L 2*° (2vr)= 



(27T) 4 * 4 U-£ ^ 



= 1 



(5.31) 



We work out this equality simplifying common factors, until we obtain an expression for 
J. We consider the barred and radiation variables to be functions of the unbarred vari- 
ables. First of all, we remind that k and k n have the same direction. We thus make the 
replacements 

d 3 k = dfl 2 k 2 dk , d 3 k n = dfl 2 kl dk n , (5.32) 

and cancel the common factor dVt 2 on both sides of eq. (5.31). Boost invariance of the 
phase-space elements and of the four-dimensional delta function also guarantees that 



rc-1 



n 



d 3 t; 



n-l 



n—1 



L 2*? (2*y 



(2tt) 4 <5 4 u-fc-j> =n 



d 3 h; 



n—1 



i=l 



i 2*? (2vr) ; 



(27T) 4 5 A [q-k n -^2k i 



i=l 



(5.33) 

In fact, k n , £, y and tp are functions of k n +i and fc only, while k±, . . . , fc n -i depend upon 
fei, . . . , k n -\ via the boost B, that depends only upon k n+ \ and k. Furthermore, according 
to (5.27), 

B (q - k) = q - k n , (5.34) 

so that eq. (5.33) is implied by boost invariance. The members of eq. (5.33) can thus be 
divided out on both sides of eq. (5.31). Performing also the replacement 



d 3 k n+ i q 2 
2k° n+1 (2vr)3 " (4vr) 



■^d^d cos tp dtp , 



(5.35) 



where ip is the angle between k n+ \ and k, on the left hand side of eq. (5.31), we can also 
divide out d£ dtp. We are left with the identity 

~ 2 k 2 dk 



(47T) 



■ £ d cos ip ' 



h0 



Jdyk n dk n , 



(5.36) 
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and, thus, we just need to express y and k n in terms of cos if) and k, at fixed £, and compute 
the Jacobian of this two- variable transformation. The following sequence of relations give 
y and k_ n as functions of cos if) and k (at fixed £) 



Kn = \lk +kn+l -2fcfc„+i COS^, 

9 



- 2 Ml, 



2q° 



Mr ec — (q° kn+l kn) i 
k 2 - h 2 - k 2 
2 k n kn+1 



(5.37) 



Thus 



dk 

dkr, 



dy_ 

dk 
dy 



where 



and we get 



so that 



d cos if} d cos if) 

k 2 = 2k n k n+1 (l-y) 



J 



kl 



rad 



ln ~2^ 

,2 x -1 



(4vr) 3 k, 



(5.38) 

(5.39) 
(5.40) 
(5.41) 



5.2.2 Inverse construction 



In this section, we describe the construction of the full (n + l)-particle phase space, given 
the barred variables 3> ra and the radiation variables £, y and <f>. We immediately have 



k° - k - F 1- 

^n+l — "Ln+1 — S 2 ' 



(5.42) 



The absolute value of the three-momentum of k n must instead be obtained by solving the 
equation for the energy conservation 



kn+kn+l + yk 2 + M 2 cc = q°, 



where 



k — \J}£i + kn+l + %kn kn+l V i 

and M 2 CC , according to eqs. (5.23), (5.27) and (5.28), is given by 

Ml c = {q-k n f. 

Under the obvious condition 



1.0 < 
K n+1 <~ 



7 2 — M 2 
1 rcc 

2q° 



(5.43) 
(5.44) 

(5.45) 
(5.46) 
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there is one and only one (positive) solution of eq. (5.43) 

q 2 - ML - 2q°k„,, 

2 [q° - kn+i (i - y)\ 

Having specified the size of k n , we construct the vectors k n and k n +i in such a way that 
their vector sum k = k n + k n+ \ is parallel to k n , and that the azimuth of fc n +i relative 
to k (the given reference direction) is <p. Having constructed k n and k n+ i, we can define 
k = k n + k n+ \ and find k TCC = q — k. We can now compute (3 according to eq. (5.26), and 
obtain 

ki = n~ 1 k i , i = l,...,n-l. (5.48) 

Thus, the inverse mapping exists provided 

q 2 - M 2 

t < , (5-49) 

and it is always unique. 

6. POWHEG in the Catani and Seymour framework 

The separation of R into singular components may require particular attention. First of 
all, if one uses the formula 

R^= Va _ R, (6.1) 

a problem may arise due to zeros in the denominator (in fact, the CS counterterms are not 
necessarily positive). We have not tried to prove that this can really happen. However, 
the problem is easily solved by using instead (for example) 



2 



R a <= \ 2 R. (6.2) 

In case the n-body cross section possesses singular regions, as discussed at the beginning of 
section 2.2, a further problem may arise. Formula (6.1) may not be adequate to partition 
the different singular components of R. This is better seen with an example. Consider 
the Z + jet production process. The n-body process corresponds to a Z + I final state, 
where I is a light parton, and the (n + l)-body process corresponds to a Z + l\ + I2 final 
state, where l\ and I2 are light partons. Consider now the counterterm corresponding 
to I2 becoming collinear to an initial-state parton. It is proportional to the underlying 
Z + li parton cross section. It is therefore also singular when l\ becomes collinear to an 
initial-state parton, since the Z + l\ cross section is singular in this limit. In standard 
NLO calculation, an infrared-safe observable O, that vanishes when two singular regions 
are approached at the same time, suppresses the singular regions of the underlying Born 
process in the counterterm (see eq. (2.19)). This problem is easily solved by writing 



fl Qr = V(«V (6 - 3) 
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where if is a positive function that vanishes when its argument approaches an n-body 
singular configuration. We again have 



and now R ai is singular only in the a r region. 

If one wishes to enforce transverse-momentum ordering (see sec. 4.1), one can first sep- 
arate R as in eqs. (2.62) and (2.63), using the S functions defined in eqs. (2.75) and (2.80), 
with the d functions defined as in eqs. (4.18) and (4.19). The contributions to R cor- 
responding to a given singular region can then be separated according to each possible 
spectator using a formula similar to eq. (6.3). 

We now discuss the various dipole configurations of the CS approach. We use the 
notation introduced in ref. [27], slightly modified in order to be consistent with the notation 
of the rest of this paper. 

6.1 Final-state singularity with final-state spectator 
6.1.1 Radiation and barred variables 

The labels k denote the radiator, the emitted, and the spectator partons respectively. 
Without loss of generality, we assume that the radiated parton j is the (n + l)th parton. 
We introduce the variable yij k, that is zero in the soft and collinear limit, 

^ ki ■ kj 2 ki ■ kj ^ 

k{ ' kj -\- [ki -\- kj^j ' k k {ki ~\~ kj -\- k k ^ 

The n barred momenta, are defined as follows 

h = z — k k , (6.6) 

1 - Vij,k 

ki = k~ij = ki + kj - — kk , (6-7) 
1 - Uij,k 

k r = k r , r = 1, . . . , n + 1, r^i,j,k. (6-8) 

Notice that we have introduced the variable kij = ki, so that kk and kij correspond to pk 
and pij of ref. [27]. The expression (6.5) for y^-^ is determined by imposing fe?- = 0, and 
since 

ki + kj + k k = k k + hj , (6.9) 

momentum conservation is enforced. 

Observe that, as required, the barred momenta satisfy the momentum conservation 
relation 

n 

k & + k e = J2h, (6.10) 

i=l 

where the initial-state barred momenta are defined as 

k® = k 9 = x^Kq (6-11) 
k e = k e = x e K e (6.12) 
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so that 

% = x© . (6.13) 

In addition, we define 

l^/Cj + Kj) ■ K k Kij ■ K k 

that, in the collinear limit, is equal to the fraction of the momentum carried by the ith 
particle. We write eq. (5.20) of ref. [27] (in d = 4) in our notation as 

(2ki ■ • k k ) dcf) 

[dki (kij,k k )] = — — — - dzi dy ijtk 6>(^(1 - 2,)) 9(y ijtk (l - y ijtk )) (1 - y^) , 

(6-15) 

where <f> is the azimuthal angle of ki in the centre-of-mass system of (kij + k k ) and 



/ 



# = 2vr . (6.16) 



The integral of a generic finite quantity F over the (n + l)-parton phase space can then be 
written as 

j dx e dx e d<& n +i (x m K & + x e -ftT e ; k±, . . . , fc n +i) £(x e , x e ) 
x (^e; ^e! ^l) ■ ■ ■ j k n +i ) 
= J ^ S dx e d^ n (x ffi -ftT ffi + x e -ftT e ; k±, . . . , fc n ) [(ifcj (fcy , £(xe,x e ) 

x ( x e > x e ! > • • • > ^n+i ) 
= j d& n d$ rad £(x e ,x e ) F(x e ,x e ;A;i, . . . ,k n+1 ) , (6.17) 

where 

(2ki ■ ■ k k ) d(j) 

d^rad = — ^2 — d *i d Vij,k 9(^(1 - Zi)) 9(y ijtk (l - y ijtk )) (1 - y ijtk ) . (6.18) 
6.1.2 Inverse construction 

In this section, we describe the construction of the full (n + l)-particle phase space, given 
the barred variables $ n and the radiation variables. The goal is to build the momenta of 
the emitted particle, kj, of the emitting particle, ki, and of the spectator particle, k k , given 
kij, k k and the three radiation variable: yij >k , z% and (j). All the other momenta remain 
unchanged, according to eq. (6.8). 

From eqs. (6.5), (6.14) and from momentum conservation (6.9) we have 

ki ■ kj = yij^ k k k ■ kij (6.19) 
ki ■ k k = Zi k k ■ kij (6.20) 
ki ■ kij = yij^ k (1 - Zi) k k • kij . (6.21) 

We make a Lorentz boost B to centre-of-mass frame of (kij + k k ), and we denote with a 
"prime" the momenta in this frame. We fix the z' axis parallel to k' k . In this reference 
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frame, the boosted momenta have the following Lorentz components 

k' k = E(1, 0,0,1) (6.22) 
% = E(1, 0,0,-1) (6.23) 
fc- = £-(1, sin 0- cos 0, sin 6»- sin 0, cos 0j) (6.24) 

where k k -kij = 2E 2 , so that 



kk ' k{j 



E = \\ J ^^ ( 6 - 25 ) 



and eqs. (6.20) and (6.21), evaluated in this frame, give 



2E 2 z { = E\ E(l- cos (6.26) 
2E 2 y ij>k (l - Zi) = E[E (l + cos 0[) . (6.27) 



Solving these equations, we derive 



E' i = \l 1 ^[y ij , k {l-~z i ) + ~z i \ (6.28) 



Uij,k{l Zi) %i 

COS (7; — — — 

yij,fe(i-^)+^ 



(6.29) 



and we have so determined all the four components of k[. We boost back to the original 
frame and we obtain 

k i = M- 1 k' i . (6.30) 
From eq. (6.6) and from momentum conservation we can then write 

h = (1 - Vij,k)h (6.31) 
kj = Vij,kkk + fcjj — ki , (6.32) 

and this completes the task of building the (n + 1) final-state momenta. 

6.2 Final-state singularity with initial-state spectator 
6.2.1 Radiation and barred variables 

In this case the radiator i is a final-state parton and the spectator k is an initial-state one, 
that we assume for definiteness to be the © parton. As before, without loss of generality, 
we assume that the radiated parton j is the (n + l)th parton. We introduce the variable 

that approaches 1 in the soft and collinear limit, and the following n barred momenta 

ki = kij = k{ + kj (1 Xij^) k^ , (6.34) 
k r = k r , r = 1, . . . , n + 1, r^i,j. (6.35) 
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Momentum conservation reads 

n 

K + k e = ^2~ki, (6.36) 



i=i 



where 



k® — x ij,®k§ — x b-K$> i (6.37) 
k e = k e = x e K e , (6.38) 



and 



Introducing the variable 



x ® ~ x ij,® x ® i (6.39) 
x e = x e . (6.40) 

(ki + kj) • k B 

that, in the collinear limit, is equal to the fraction of the momentum carried by the ith 
particle, we can write the integral of a generic finite quantity F over the (n + l)-parton 
phase space as 



dx B dx e d<§> nJr i (xqKq + XqKq] k±, ... , k n +i) £(x ffi , x e ) 
x F { X (B i x e i k\ i ■ ■ ■ , fen+l ) 

x F(x 9 ,x e ;ki,...,kn +1 ) , (6.42) 
where [dki(kij] k e ,x)] is given by eq. (5.48) of ref. [27] in d = 4 dimensions 

[dki{kij;k B ,x)] = ^ 2k ^ ^ dzi dx ij:@ 0(^(1 - Zj)) 6(x(l - x)) 5(x - x ij:@ ) , (6.43) 

and 4> is the azimuthal angle of ki in the centre-of-mass system of (kij + k B ). 

Performing the integration in x, that fixes with the change of variable of 

eqs. (6.39) and (6.40), we obtain 

/dx ■ ■ _ _ fx \ 

cte e dx e — ^ d$ n (s e ^ + x e K e ; ki . . . k n ) C I — x e ) 

(2k l3 ■ k B ) d<P_ ^ (l _ ^ _ ^ 



X - 


167T 2 


x e 


{ x ij,@ ~ 




( x ® 


'rad £ 


( ,i 







,x e F(x e ,x e ;A;i,...,A; n+ i) , (6.44) 



where 



d *«d = (2 ^ ffi) ^ ^ C 1 " *)) " xy, e )) ^ (* - x e ) . (6.45) 
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6.2.2 Inverse construction 

In order to reconstruct the (n + 1) momenta from $ n and $ ra d; we follow closely the 
procedure used in section 6.1.2. More precisely, we have to reconstruct ki, kj and k 9 given 
k^, k B and the radiation variables Xy )ffi , Zi and <p. 
From eq. (6.37) we immediately have 

K = , (6.46) 

and using eqs. (6.33) and (6.41) we can write 

ki ■ kj = (1 — Xij t(B ) • kij (6.47) 

ki ' = ^i fee • kij (6.48) 

ki ■ k^ = (1 - x ijt@ ) (1 - Zi) k k ■ kij . (6.49) 

These equations are similar to eqs. (6.19)-(6.21) with the substitutions 

k k k 9 (6.50) 

Vij,k <-> 1 - x ijiB (6.51) 

so that, in the centre-of-mass of the (kij + fc e ), the energy and the angle 9[ of ki with the 
k' direction are given by (see eqs. (6.28) and (6.29)) 



El = yj [(1 - xy, e )(l " Si) + H (6-52) 

cos 6'i = = (l-^e)(l-^)-^ (6 _ 53) 
(l-^, e )(l-^) + ^ 

We can boost back to the original frame and obtain ki. The last momentum kj is con- 
strained by momentum conservation 

kj = kij — ki + (1 — Xjj ie ) /c ffi . (6.54) 

6.3 Initial-state singularity with final-state spectator 
6.3.1 Radiation and barred variables 

In the case where the parton i (assumed here to be the (n + l)th parton) is emitted from 
an initial-state parton, that we take for definiteness to be the © one, and in the presence 
of a final-state spectator k, we define the following n barred final-state momenta 

kk = h + k k - (1 - x ikySj ) k m (6.55) 
k r = k r r = 1, . . . , n + 1, r / i,k , (6.56) 

where the variable ^ ^ 

(ki + k k ) ■ k m 
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approaches 1 when fcj becomes soft, and approaches z in the collinear limit, when k{ = 
(1 — z) k @ . Momentum conservation reads 

n 



i=i 



where 



~~ x ik,® ~~ X (B-K® (6.59) 
k e = k e = x e K e (6.60) 



so that 



Introducing the variable 



X S> ~ x ik,® X (B (6.61) 

x e = x e . (6.62) 

•'-(*. <6 - 63) 

we can write the integral of a generic finite quantity F over the (n + l)-parton phase space 
as 

J dx @ dx e d® n+1 (x ffi i^T ffi +x e K e ;ki, . . .,k n+1 ) £(x ffi ,x e ) 
x F ( x ®i x e] k\, . . . , k n +i) 

x F(x ffi ,x e ;fci,...,fc n+ i) , (6.64) 
where [<i/cj(/cfc; fc e , x)] is given by eq. (5.72) of ref. [27] in d = 4 dimensions 

[dki(k k ;k s ,x)] = ^ dx ik ^ 9(ui(l - u,)) 6>(x(l - x)) <5(x - x ifcie ) , (6.65) 

and </> is the azimuthal angle of k% in the centre-of-mass system of (k k + k B ). 

Performing the integration in x, that fixes x = Xi kj(B , with the change of variable of 
eqs. (6.61) and (6.62), we obtain 

j dx e dx e dXtk >® d$ n (x @ K 9 + x e K e ;ki ...k n )C ( -^-,x e ) 

(2k k • fc ffi ) # f . .x / , 

— i&Tr^ — 2vr ^V u n ~ ^(^M 1 ~ x ik,®)) 

x (x«fc,e ~~ ^e) F (x ffi , x e ; k±, . . . , fc n +i) 
= / d* n d$ rad £ ( -^-,z e ) F(x ffi ,x e ;/ci,...,A; n+ i) , (6.66) 

J \ x ik,e / 

where 

d$rad = 2^ ^T" dUi d ( U i( 1 ~ U i)) d { x ik,®{ 1 ~ x ik,(s)) # (x ifcjffi - X ffi ) . (6.67) 

107T Z7T x ik,<$ 
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6.3.2 Inverse construction 

This case is completely analogous to the one in section 6.2.2 with the substitutions 

kj <-> kk (6.68) 

kij <-> k~k (6.69) 

x ij,® x ik,® (6.70) 

Zi <-> Ui (6-71) 



so that 



fc® = (6.72) 



E[ = \j [(1 - x lK(B )(l - ui) + m] (6.73) 
co B (? i = = ?- X *>°l?- Ui) - Ui . (6.74) 

where E[ and 9[ are the energy and the angle that the vector k\ forms with the direction 
of k' 9 , in the centre-of-mass of the (kk + A; e ). The four- vector fc, is obtained from k\ with 
a boost back in the original frame, while kk is constrained by momentum conservation 

kk = — ki + (1 — Xik,®) • (6.75) 

6.4 Initial-state singularity with initial-state spectator 
6.4.1 Radiation and barred variables 

In the case where the parton i (that we take to be the (n + l)th parton) is emitted from an 
initial-state parton, that we take for definiteness to be the © one, and in the presence of an 
initial-state spectator, the parton, we define the following n barred final-state momenta 

k£ = A>*„(K,K)k» r = l,...,n + l, r^i, (6.76) 

where the boost tensor is given by 

A» V {K, K) - g% {K + ICy K 2 ' ( ' 

n+1 

K = k S) + k e -k i = ^2k r , (6.78) 

r=l 
n+l 

K = Xi t9e k @ + k e = k r , (6.79) 



r=l 



and 



fe ffi • k e 
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Notice that, in the soft limit, Xj iffie approaches 1, while, in the collinear limit, i.e. fej 
(1 - z) fe e , it approaches z and that F> = A^ U (K, K) K u , so that K 2 = K 2 . 
Momentum conservation reads 



fee + h = ~ ki ( 6 - 81 ) 



i=i 



where 



so that 



fee ~~ x i,®e ^© ~~ X ®K® (6.82) 
fee = fee = x e K e (6.83) 



x ® — x i,®e x ® (6.84) 
x e = x e . (6.85) 



Introducing the variable 
we can write the integral of a generic finite quantity F over the (n + l)-parton phase space 



as 



j dx m dx e d<& n +i (x§K 9 + x e K e ; k±, . . . , fe n +i) ^(%! ^e) 

x F ( x ® > x e i > • • • ) fera+l ) 
= / ^ cfc ». (,^. + * e * e ; . . . I,) [* (*., * e , x,] , e ) 

x F (x ffi , x e ; fei, . . . , k n+ i) , (6.87) 

where [dfej (fe e ,fe e ,x)] is given by eq. (5.151) of ref. [27] in d = 4 dimensions 

[dfei(fe®, fe e , x)] = ^^a 9 ^ ^~ ^ dx i,®e e 61 ^ ~ ~ ^ ^ ~ Xi '® 9 ^ ' 

(6.88) 

and 4> is the azimuthal angle of fej in the centre-of-mass system of (fe e + fe e ) or in the 
laboratory frame. 

Performing the integration in x, that fixes x = Xj, ffi e> with the change of variable of 
eqs. (6.84) and (6.85), we obtain 

/dx ' — — ( x \ 

dx 9 dx e — d$ n (x e K e + x e K e ; ki . . . k n ) £ I — — , x e ) 

m ' ! ] # dvi e e (i - r -^— (?(x i>ee (i - x i>ee )) 

V ± x i,<se J 



16vr 2 2vr 



= / d$ rad £ ( ^^,X e J F(x ffi ,X e ;fel,...,fen+l) , 



(6.89) 



- 55 - 



where 

d*»d = g dv t ^ (*) (l - T -|— ) *(s i>ee (l - x i>ee )) 

x6>(x iiffie -x ffi ) . (6.90) 

6.4.2 Inverse construction 

To reconstruct the full (n + l)-particle final state, given fc e , fe e (i.e. x ffi and x e ), the n 
final-state momenta k r and the three radiation variables Xj i9e , £j and 0, we follow the 
same procedure already used in section 6.1.2. From eqs. (6.82) and (6.83) we immediately 
have 

^ = (6-91) 
k e = k e (6.92) 



that is 



We compute then 



Xi 



x e = x e . (6.93) 



ki- k @ = Vi k 9 ■ k e (6.94) 
ki • k e = (1 — Xi j(f) Q — Vi) • /c e (6.95) 

in the centre-of-mass system of (/c e + k e ), and we find that the energy E[ and angle 6[ that 
the vector k\ forms with the k' direction are given by 



E[ = f-^{l-x h9e ) (6.96) 

1 ^1,06 

We boost k[ back in the laboratory frame and we obtain ki. Once ki is known, we can 
build the two vectors 

K = k m + k e -ki , (6.98) 

K = k B + k e , (6.99) 

and the inverse of the boost tensor A 

a^,^- ^^ ^, (6 ,oo) 

and compute the remaining n momenta 

k r = A-\K,K)k r , r = l,...,n. (6.101) 
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7. Examples 

In this section we show in detail how to implement the POWHEG formalism in two simple 
cases: e + e~ — > qq and vector-boson production at a hadron-hadron collider, i.e. h m h e — > V. 

7.1 e+e _ — > in the Catani and Seymour formalism 
Born, virtual and real corrections 

We consider the scattering e~{p\)e + {p2) — > 7* — ► q{k\) q(k2) g{k 3 ) where, in the centre- 
of-mass frame, 



pi = y(l,0,0,l) (7.1) 

p 2 = ^(1,0,0,-1) (7.2) 

fci = fcj (1, sin 61 cos 0i, sin 9\ sin 0i, cos #1) (73) 

&2 = ^2 iXi s i n ^2 cos 02, sin 62 sin 02, cos 62) , (74) 

and define 

2 (7 • /cm 2 3} 

q = Pi+P2, k 3 = q-k 1 -k 2 , £{1,2,3} = 2 ' ' • ( - 7 ' 5 ) 

The Born cross section 5, in arbitrary units, is given by 

B (h, h) = B (h) = 1 + cos 2 h , (7.6) 

where k\ and k,2 are the momenta of the quark and antiquark at Born level, and are given 
by 

k\ = — (l, sin 61 cos (pi, sin 9\ sin <fii, cos #1) (7.7) 

/c2 = — (l, — sin 6*i cos <pi, — sin6>i sin 0i, — cos #1) . (7.8) 

In this section we drop the flavour index in the Born cross section, for ease of notation. 
In virtue of the azimuthal symmetry of the cross section, we could fix from now on (f>\ = 0, 
and perform a random azimuthal rotation at the end of the generation. Here we keep it 
variable. The Born two-body phase-space element is given by 

so that $2 = {0i,4>i}, while the finite soft-virtual contribution of eq. (2.107) is equal to 

V(e 1 ) = -C F B(9 1 ). (710) 

The real-emission term R 

_ 8vrC F a s x\ (l + cos 2 6>i) + x\ (l + cos 2 2 ) n , , 

-R= 2 T\ ui ^ l 7 - 11 ) 

r (1 - ari)(l - x 2 ) 

has two soft collinear regions. We call region 1 the k\-k% — > (x 2 — > 1) region, and region 
2 the fc 2 • /C3 — ► (xi — ► 1) region. 
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Counterterms 

The counterterms can be computed using eqs. (5.2) and (5.7) of ref. [27] and are given by: 
1. Region 1 (ki ■ k 3 0): = {y 13 ,2, z v 0} 



cw^ 1S2= 87rC r s — 



¥ 2/13,2 



l - 5i (l - 2/13,2) 



where 



&1 = &13 = 9 



"&2 , 



fc 2 



B(fci,fc 2 ) (7.12) 
-*2, (7.13) 



1 - 2/13,2 1 - 2/13,2 

so that 

01 = 7T — #2 , 01 = (tT + 02 ) mod 27T , 

since fc 2 is only rescaled with respect to fc 2 - 

We can express the counterterm in terms of the x\ and x 2 variables too, and we 



(7.14) 



have 



C (i) 



87rC F a s 



1 — x 2 V 2 - xi - x 2 



(1 + si) + 



where we have used eqs. (6.5) and (6.14) 
2/13,2 = 1 - x 2 , 

so that 

h = q , 

X2 



1 — Xl 

X2 



_ Xl + X 2 - 1 

X2 

k 2 = — ■ 

%2 



B(k u h), (7.15) 



(7.16) 



(7.17) 



This way, we can use two sets of radiation variables for the first singular region: the 
{2/13,2, z v 4>} set or {xi, x 2 , 4>] one. 



2. Region 2 (k 2 • k 3 0): = {2/23,1, z 2 , 0} 



C ( 2 )^P 23il = ^J_ 



where 



fci 



9 2/23,1 

-hi , 



1 — 5 2 (1 — 2/23,1) 

^2 = k 2 3 = q 



(1 + S2) 



1 - 2/23,1 " 1 - 2/23,1 

so that 

0~i = 6>i , <j>i = 4>i, 

since fei is only rescaled with respect to fei. 
In terms of the xi and X2 variables, we have 



ki , 



(7.18) 
(7.19) 
(7.20) 



8irC F a s 



1 

1 — Xl V 2 — Xl — X2 



-(l + x 2 ) + 



1 — x 2 

Xl 



Bfiuh) (7.21) 
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where we have used eqs. (6.5) and (6.14) 

1/23,1 = 1-*!. z 2 = Xl+X2 ~ 1 , (7.22) 

X\ 

so that 

ki = ~, h = h3 = q--. (7.23) 

X\ X\ 

This way, we can use two sets of radiation variables for the second singular region: 
the {2/23,1, z 2 , 4>) set or {xi, x 2 , 4>] one. 

Radiation phase space and inverse construction 

The radiation phase-space element of eq. (6.18) has the same functional form in both the 
two regions: 

^rad = dz dy (1 - y) 6(z(l - z)) 6(y(l - y)) (7.24) 

dxi dx 2 0(1 - xi) 9(1 - x 2 ) 6(xi +x 2 -l), (7.25) 



16vr 2 2vr 

and we can use either of the two sets of radiation variables: {y, z, <ji\ or {x\,x 2 , </>}, where 
we identify y = 1/13,2, z = z\ in region 1, or y = 1/23,1, z = z 2 in region 2. In each of the two 
regions, we can express the angles 0\ and 9 2 in the real term (7.11) in terms of the Born 
barred variables, Q\ and (fi\, and in terms of the radiation variables, that here we take to 
be x±, x 2 and (j), once all the 3 final-state momenta have been constructed, according to 
the procedure described in section 6.1.2. For example, considering the first singular region 
and the two Born variables 6\ and 4>i (that fix uniquely k\ and k 2 ), we have 

2 = vr - 0i , 4> 2 = (fa + tt) mod 2vr , (7.26) 

and k[ of eq. (6.24) is characterized by 



E[ = ^f Xl , (7.27) 

C os9[ = 2 - 2Xl - 2X2 + XlX \ (7.28) 
X\ x 2 

and by a random angle (j). Observe that, in this case, no boosts are needed for the inverse 
construction, since the dipole rest frame coincides with the e + e~ centre-of-mass frame. 
The vector k\ is simply obtained as follows 

ki = Rz(h)Ry(S 2 )K, (7-29) 

where R z / y (a) is a rotation around the z/y direction of an angle a. The remaining two 
momenta are then given by 

k 2 = x 2 k 2 , fc 3 = q — ki — k 2 . (7.30) 

The corresponding results for the second region are obtained from the previous ones with 
the exchange of indexes 1 <-> 2. 
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Generation of the Born variables 

The B function of eq. (4.13) is then given by 

\rM _ c(i) 



B{e l )=B{e l )+v{h) + / d$ rad 



R m _ C (2) 



where 



RW = R- 



(7.31) 



(7.32) 



CW + ' + C( 2 ) ' 

In order to generate the two Born variables #i and </>i distributed according to eq. (7.31), 
we need to construct the B function of eq. (4.23). We then make the following change of 
variables 



X2 



2vrX 



(3) 
rad 



SO that 



and we obtain 



d<f> rad = 



Q X {1) dX {l) HX {2) JX {3) 

J£~2 A rad « A rad « A rad « A rad 



(7.33) 

(7.34) 



B(S 1 ,x nd ) = B(5 1 ) + v(e 1 ) + 

that satisfies 



16vr 2 



X 



(i) 

rad 



(rW-CW) + (rW-CW)], (7.35) 



fdx^ fdx® f 1 dx^ d B(e 1 ,x rad )=B(e 1 ) 

Jo Jo Jo 



(7.36) 



The function B can now be integrated over the full 3-body phase space, using an integrator 
that can generate unweighted events, like the SPRING-BASES package, so that one can 
generate the Born configurations very efficiently. 

Generation of the radiation variables 

Once we have generated the Born configuration distributed according to the B function, 
one must generate the hardest radiation. We first define the k T of the radiation. The 
definition is ambiguous to some extent. The only requirements that it must satisfy is that 
it must be of the order of the radiation transverse momentum in the collinear limit, and 
that it must coincide with it in the soft-collinear limit. A suitable definition is 



k% = q 2 yx 3 = q 2 y [1 - z(l - y)\ , 

where X3 is the energy fraction of the gluon. 

The Sudakov form factor of eq. (4.16) is given by 



(7.37) 



A($ 2 ,Pt) = exp< 



exp 




B{K) ^t-PtJ 



(7.38) 
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In order to apply the veto method of appendix A to our case, following the procedure of 
section 4.3.2, we need to specify the form of an upper bounding function. A suitable choice 
turns out to be 



Na s (kl) 



1-2/ 



y l-z(l-y) 



E 



(1-9) 



where 



Qre{i,2} 
32ttC f 



N 



(7.39) 



(7.40) 



We denote the integral of the upper bounding function by 

A u (p T )=exp{-^ Pg. Cdz Cdya s {kl) l^L- ± Mk^-pl) 

\ it J 2vr J Jo y 1 - z(l - y) 

The integral in the exponent can now be easily performed 



A„(p T ) = exp { - 



2C, 



7T Jo Jo k 2 /q 2 
and introducing the dimensionless quantities 

~k 2 T = ^=y[l-z(l-y)}, 

and the integration over dk 2 , we have 



n 2 

~2 _ Pt 

Pt 9 1 

q 2 



(7.41) 
(7.42) 

(7.43) 



A u (p T ) = exp | dk% jf 1 dz £ dy (1 - y) 5 (y (1 - z(l - y)) - fc 2 ) 



(fc 2 g 2 ) ] 



1.2 



exp < 



exp 



dk 



2C P r 1 70 a A k W 



71 H\ 



2C F Z* 1 ~o as \ k ' q 



^fdk 

11 JP 2 T 



dy 

/,■? .//.••: y 

2J2 1 



U2 



1 



The integral can now be easily performed using the one-loop expression 

1 



a s {k 2 T ) 



6 log(fc2/A2) ' 



and we obtain 



A u (p T ) = exp 



p 2 g 2 log ( g y A 2 ) 

l0g 7 + l0g A^ l0g l0g(^) 



(7.44) 



(7.45) 



(7.46) 



Observe that we have a lower limit on the acceptable values of p T , since a s (p T ) must be 
well defined for both the two- loop and the one-loop expression of a s . We thus introduce a 
P™ m ~ A. Summarizing: 
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1. Set p 



max 
T 



q 2 , where q 2 is the maximum value of p 2 , such that A u (p 



max \ 
T J 



1. 



2. Generate r with < r < 1, and solve the equation r = A u (p T )/A u (p™ ax ) for p T . If no 
solution is found, or if p T < p" im no radiation is generated, and the event is returned 
as a qq event. 

3. If a solution with p T > p™ in is found, the radiation variables, according to the veto 
technique described in appendix A, are then distributed according to 



1-1/ 



1 



y l-z(l-y) 



% 2 y(l-z(l-y))-p2). 



(7.47) 



The <5 function poses no constraints on the azimuthal angle 4>, that is then generated 
uniformly between and 2ir, but it fixes z to be 



v - pl/q 2 
y(i-y) 



z = 



(7.48) 



where the probability distribution of y is dy/y, so that y is uniformly distributed in 
logy, between the range imposed by < z < 1, 



log ('^) < logy <^ log 



4. Generate a uniformly-distributed random number r 1 in the range 

1-y 1 



(7.49) 



< / < Na s 



If 



r'< 



E 

a r €{l,2} 



y l-z(l-y) 
R{6i,4>i,4>,z,y) 

B{h) 



(7.50) 



(7.51) 



accept the event. Otherwise set p^ ax = p T and go to step 2. 
5. If 



r' > 



(i-y) 



(7.52) 



then set a r = 2. Otherwise set a r = 1. 
This completes the generation of the radiation variables and of the singular region a r . 

7.2 e+e - — ► qq in the Frixione, Kunszt and Signer formalism 
Virtual and real corrections 

Using the notation of section 7.1, we write the virtual corrections, normalized as in eq. (2.92), 
as 



V b = N 



2vr 



C F 3C F \ 2C F ^ 

"2 I — + ) H log — 7T + V fin 



2e 



B{h) , 



(7.53) 
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where 



V fin = C F - 8 + 3 log ^ " log 2 ^ J (7.54) 
and 

B{h) = l + cos 2 #i . (7.55) 
Setting £ c = 1 and 5 Q = 2, equations (2.100) and (2.101) become 



13 2vr 2 \ 3 , q 2 



(7.56) 



Q = 2C F 

r u = ilog»£-£, (7.57) 



and, using eqs. (2.97) and (2.99) we get 



V(h) = J [S + 2 C F T 12 + V fin ] , (7.58) 



where the dependence upon the scale Q completely cancels. We now consider the real 
emission term R. Following section 2.4.1, we define the functions dij = (k® kj) a (1 — cos 9ij) b 
(see eq. (2.68)), so that 

d 31 = 2 b - 2a {q 2 ) a (x lX3 ) a - b (l - x 2 ) b , (7.59) 
d 32 = 2 b ~ 2a (q 2 ) a (x 2 x 3 ) a - b (l - Xl ) b , (7.60) 

where a and b are positive real numbers. We introduce the functions <Sy (see eq. (2.75)) 

l 

Sij= -rfx - ue {31,32}. (7.61) 

(I31 d 32 

Observe that there is no need to include the h factor of eq. (2.75), since only parton 3 has 
an associated soft singularity, and therefore we only consider the regions 31 and 32 with 
h = 1 (i.e. we do not consider the regions 13 and 23). Furthermore, the region 12 does not 
correspond to any singularity of R, and thus it is not included. 
We define (see eq. (2.81)) 

c 2k 3 - fo-h k 3 - k 2 . . 



7,0 7,0 ' M 7.0 7,0 ' 

"-3 "a n-3 ^2 



and (see (2.85)) 



R = R 31 + i? 32 , R 3{12} = | (|) ( l- J 3{12} ) " ^ ' (7 - 63) 



+ 



where 



R ij = RS ij , ije {31,32}. (7.64) 



The underlying Born configuration for each singular region is easily identified as the one 
that preserve the direction of the non-collinear parton. Thus, the underlying Born config- 
uration for the 32 region has ki oc ki, and for the 31 region has k 2 oc k 2 . 
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Radiation phase space and inverse construction 

Using eqs. (5.31) and (5.41), we can write the three-body phase space in terms of the barred 
and radiation variables. We work out explicitly the formulae for the 32 singular region. 
The 31 one is treated in the same way. We have 

d$ 3 = J 32 d£ dy 32 d(f>d$2, (7.65) 

where (see eq. (7.9)) 

d§ 2 = —L .(fcostfi d4>i . (7.66) 

In virtue of the azimuthal symmetry of the cross section, we can fix from now on <fii = 0, 
remembering that, at the end of the procedure, we should rotate the whole event by a 
random azimuthal angle. From eq. (5.47) we get immediately 

2(1-0 , (i + feK nKn 

and, using the fact that x\ + x 2 + £ = 2, we also get 

-, (1-1/32)^(1-0 ,- Ra x 

X1 = 1 ~ 2-ed- y32 ) • (7 ' 68) 

Given 9±, £, 7/32 and <f> we can fully reconstruct fci, k 2 and £3. Finally, from eq. (5.40), we 
get 

Generation of the Born variables 

The B function of eq. (4.13) is given by 



B(e 1 )=B(9 1 )+V(9 1 ) + J2 J WdtdyJR(h,<i>,Z,v) 



a r G{32,31} 



(7.70) 



where we are using the "context" notation introduced in eq. (2.21). We now explicitly 
show how to deal with the distributions, by illustrating the computation of the 32 term as 
an example: 

J di dy 32 d<P J 32 i?32 = 7^3 J d£, dy 32 d<t> (1 " m2) ? ^ 32 l^f ' 

(7.71) 

One can easily verify that both J32 and (1 — 2/32) £ 2 -R32 have a finite limit for 2/32 1 or 
£ — > 0. We can thus expand the distributions according to their definitions (2.88), (2.89) 
and (2.90). Introducing, for ease of notation, the function 

a 32 (01, <i>, e, 2/32) = (1 - ys2) e R32 , (7.72) 
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we have 



2tt 



d<p 



f d£ [ dy 32 J32 R32 = jj-yT [ d(f> [ dt, [ dy 32 7 — — — 

J-l ( 47r ) Jo JO J~l i 1 - 2/32 



^32(ei,^C,y32)-(l-0 M2{hA,i, 1) 



452 (01, </>, 0,1/32) -^32(^1,^0,1) 



(7.73) 



We observe that the Jacobian of eq. (7.69) has an integrable divergence at (£, 1/23) — ► (1, —1) 
region. This happens because the recoiling system is massless, so that there is no upper 
bound on £ (see eq. (5.49)). In order to clarify the nature of this singularity, we introduce 
the parametrization 

£ = l-pcos77, y 2 3 = -1 + psinr], (7.74) 



and from eq. (7.69) we get 



J: 



32 



4(1 -m 



(47r) 3 (2(l-£)+£(l + 2/32)) 2 (4vr) 3 p 



4 cos 7/ 



(2 cos r\ + sin 7])' 



+ o( P ) 



(7.75) 



diverging as 1/p in the small p limit. This singularity is integrable, but cannot be inte- 
grated using Monte Carlo techniques, since it yields a divergent error. On the other hand, 
<^32\/l ~~ + 2/32 has no divergence, which suggests the use of square-root importance 
sampling for both 1 — £ and l + t/32 in order to overcome the problem. We thus parametrize 
the integration variables as 



« = i-(-^) 2 . 



y23 



-1 + 2 



( X rad) ' 



(7.76) 



The generation of the Born variables is performed along the lines of section 4.3.1. We 
define 



a r e{32,31} 



so that 



Jo Jo 



dX, 



(2) 







(7.77) 



(7.78) 



Then one performs a full integration of B, using an integrator that can generate unweighted 
events, like the SPRING-BASES package. After this step, one can generate Born configura- 
tions (i.e. 9{) distributed according to B. 

Generation of the radiation variables 

After having generated the Born configuration distributed according to the B function, we 
must generate the hardest radiation. We first define the k T of the radiation. As we have 
already pointed out, the definition is ambiguous to some extent. The only requirements 
that it must satisfy is that it must be of the order of the radiation transverse momentum in 
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the collinear limit, and that it must coincide with it in the soft-collinear limit. A suitable 
definition is 

2 



(7.79) 



The no-radiation probability, according to eq. (4.16), is given by 



A(0i,p T ) = exp 



k a r e{32,31} 



/ 



d(f) d£ dy- 



B{h) 



6{k T — p T ) 



(7.80) 



In order to generate the radiation variables distributed according to eq. (7.80), we make 
use of the veto method of appendix A. We need then to specify the form of an upper 
bounding function. A suitable choice turns out to be 



Na s (k T ) 



E 



arG{32,31} 



B(B 1 ) 



(7.81) 



where y stands for the common value of y%i and 2/32 on the right hand side. Notice that 
U{^,y) is singular also for y — > —1. In this the upper bound properly covers also the 
singularity in J for y — > —1, (see eq. (7.75)). The normalization iV is evaluated by probing 
the 9i,(f>, £, y phase space randomly with an adequate number of points. The a s (k T ) factor 
appearing in eq. (7.81) corresponds to the one-loop expression. 18 
We now want to generate (ft, £, y according to the distribution 



exp 



dtfdt'dy' u(e,y')e(k' T -k T ) 



U&y) d<\>didy, 



(7.82) 



where k' T = k T (£',y') and k T = k T (£,y), and then use the veto technique of appendix A 
to correct to the exact distribution. One generates first the k T value. Its probability 
distribution is given by 



exp 



b'di'dy'U{i\y')e{k' T -k T ) 



U(£,y) dcj)d£dy5(k T - p T ) 



dA u (p T ) 



dpi 



where 



exp 



d(/ ) d^dyU^,y)9(k T -p T ) 



(7.83) 
(7.84) 



The integral in the exponent of eq. (7.84) is easily manipulated to yield 



dcj) d£ dy 



Na s (k T 
(l-2/ 2 U 



lax 

6(k T — p T ) = 2ttN I — - a s (k T ) log 



PT 



1 + yi - (fc T /fc max ) 2 

1 - y/l - (k T /k n 



\2' 



with 



g u /2. 



(7.85) 



(7.86) 



Since the e + e — ► qq cross section is of zeroth order in the strong coupling constant, the use of one-loop 
as in the radiative corrections is adequate in this case also in the full NLO formula. 
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So, we want to generate p T according to the Sudakov form factor 



exp 



-2vriV 



rk-a 
JO 



dk T 1 + y 1 {k T /k mSLX ) 2 

— a s {k T ) log- v ==g(fc T -p T ) 

fc T 1 - Vl - (fcxAmax) 2 



(7.87) 



which has again the form of eq. (7.84) for a single variable k T instead of the set of variables 
£, y. The integral in eq. (7.87) is still too complex to be performed analytically, so we resort 
a second time to the veto method. Using the inequality 



W > j 1 + VI ~ (fc T /fc ma ^ 

" 1 - VI - (^max) 2 ' 
we generate the p T distribution according to dA uu (p T ), where 

f^max A I. A / , 2 



(7- 



A„w(p T ) = exp 



-2ttN 



exp 




^a s (fc x ) log%^(fc T -p T ) 



4k 2 



log !2il^)_ log .^' 



log(p2/A2) 



(7.89) 



Observe that we have a lower limit on the acceptable values of p T , since a s (p T ) must be 
well defined for both the two- loop and the one-loop expression of a s . We thus introduce a 
Pt m ~ A- Summarizing 



1. Set p; 
A un (p; 



max 
T 

maxj 



k n 
-- 1. 



the maximum allowed value according to eq. (7.86). We have 



2. Generate a uniformly-distributed random number r between and 1, and solve the 
equation r = A MU (p T )/A MU (p™ ax ) for p T . If no solution is found, or if p T < pi 
there is no radiation, and a qq event is generated. 



mm 
T ) 



/ / Ak^ 

3. If a solution with p T > p" im is found, generate r with ^ r ^ log — J § a2L . 



If 



r' < lot 



max f 



go to step 4. Otherwise set p, 



max 
T 



p T and go to step 2. 



1- V /l_ (PT/femax) 2 

4. At this point p T is generated according to eq. (7.87). The radiation variables </>,y,£ 



are then distributed with a probability proportional to 

1 



:<*(M£>2/) - Pt) d<j) dy d£ , 



(7.90) 



£(i-y 2 ) 

where fc T (£, 2/) is given by eq. (7.79). The 5 function does not constraint the azimuthal 
angle 4>, that is then generated uniformly between and 2tt, but it fixes £ to be 

2p T 



^ (? o v /7^2' 

and, upon integrating in £, y is distributed with a probability proportional to 



(7.91) 



V q J 1 ~y 



(7.92) 
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One thus generates a uniform random number r y between the minimum and maxi- 
mum value of the logarithm allowed by the function 



log 



1 + \/l - (Pr/fcmax) 2 " 
1 - \A ~ (PrAmax) 2 



< r„ < lot 



1 + \/l - (PrAmax) 2 " 
1 - \A ~ (PrAmax) 2 



and solves the equation 



for y. 



log 



i + y 



(7.93) 



(7.94) 



5. Now the last veto: generate a random number r" between and U(£,y). If 



a r G{32,31} 



accept the event. Otherwise set p 
6. If 

r" > 



max 
T 



B(h) 
p T and go to step 2. 



(7.95) 



jR{hA,Z,v) 



B(h) 



a r =32 



set a r = 31. Otherwise set a r = 32. 



This completes the generation of the radiation variables and of the singular region a r . 
Observe that, in our simple case, there is total symmetry between the regions 31 and 32. 
One could simply double the contribution of region 32, and, at the end, choose one of the 
two regions with equal probability. The rather pedantic discussion given above has the 
only purpose of better clarifying the method in general. 

7.3 h & h e — > V in the Catani and Seymour formalism 
Born term 

In this section, we consider the production of a massive vector boson V of mass M in 
hadron collision, in the CS formalism, 



V(k 



(7.96) 



with k\ = M 2 . Denoting, as usual, with S the centre-of-mass energy of the two incoming 
hadronic beams, S = (K @ + K e ) 2 , and with Y the rapidity of the produced vector boson, 



- 1 fc° + l\ 
2 S P — P ' 



(7.97) 



we have, from eq. (2.7), 



d 3 h 



2vr 



(27r) 3 2fc? 



(2tt) 4 5 4 (k @ + k e -k 1 ) = —dY, 



(7.98) 
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with the constraints 



and momentum conservation implies 

2 fc ffi • fe e = x^XqS = M 2 . (7.100) 

The only Born variable is then Y, so that 

*! = {Y} . (7.101) 

The squared tree-level amplitude, averaged over color and polarization of incoming partons 
(l/(4iV c 2 )), is given by 

\M B {k 9 ~ke)\ 2 = ^9 2 (*e ■ *e) , (7.102) 

where g is the Vqq coupling. Introducing the flux factor 1/(2M 2 ), we obtain the differential 
cross section 

B M -(*i) = B qq (k @ ,k e ) =^f c 9 2 , (7-103) 

independent from the external momenta. We denote by qq the Born flavour index /;,, that, 
in this example, carries information on the flavour of the (massless) incoming partons along 
the © and e direction. Using a standard convention for the numbering scheme, we have 
q = —5, . . . , —1, 1, . . . , 5, and q = —q. 

Virtual corrections 

With the scale choice Q 2 = M 2 , we have from eqs. (2.92) and (2.107) 

V fi „, 99 ~(*i) = C F (vr 2 - 8) B M -(*i) , V M -(*i) = ^ Cp M*0 • (7-104) 

Real corrections 

Three different sub-processes contribute to the real radiation production 

q(k 9 )q(k e ) - m)5(fc 2 ), (7.105) 
q(k @ )g(k e ) - V(fci) <?(**) > (7.106) 
g(k @ )q(k e ) - V(fci)9(&2)- (7.107) 

Introducing the usual Mandelstam variables 

s = (A; ffi + fc e ) 2 = (fci + fc 2 ) 2 , (7.108) 
t = (fe ffi - fci) 2 = (fc e " k 2 ) 2 , (7.109) 
u = (^ " k 2 f = (k e - h) 2 , (7.110) 

we can express them in terms of the two sets of radiation variables (see section 6.4): the 
© set, <3?® ad = {x me ,v m ,(f>}, that describes the emission from the © parton, and the © set, 
^rad = { x ee> %) 0}) that describes the emission from the © parton 
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rad 



M 2 



t — s (1 x ee u e ) 



w = — s v ff 



(7.111) 



rad 



M 2 



t = -sv e , 



u= -s (1 - x ee - v e ) , 



(7.112) 



where we have dropped the i subscript, since we have only one radiated parton, and we 
have put a © subscript on the angular v variable, in order to distinguish the two sets. 
The squared amplitudes corresponding to the processes in eqs. (7.105)-(7.107) are: 

1. The squared amplitude for the q(k m ) q(k e ) — > V(ki) gfa) process, averaged over 
color and polarization of incoming partons (1/(4N^)), is given by 



|-^(J(?| 2 — 2 ^ x y s 



F 2 2 
5 5, 



u t 2sM 
t u tu 



21 



o ^ F 2 2 

2 — q qt 



v 



® 



1 x me v@ 



+ 



1 ^mp> Va 



2xa 



(7.113) 
• (7-114) 

This expression has singularities when the gluon is emitted both along the ® and the 
e direction {v @ — > 0). 

2. The squared amplitude for the q(k m ) g(k e ) — > V(fei)(/(fe2) process, averaged over 
color and polarization of incoming partons (1/(4N C (N 2 — 1))), is equal to 



\M gg \ 2 



Tf 



-2^9 2 9 2 s 



s t 2uM 
t s st 



21 



2 |f 5 2 5s — [v% + 1 - 2x ffie (1 - x ee - e )] . 



(7.115) 
(7.116) 



3. The squared amplitude for the g(k @ ) q{k e ) — > V(A;i)(j(fe2) process, averaged over 
color and polarization of incoming partons, is given by 



Tf 

N r . 



-2^9 2 9 2 s 



u s 2tM 
- + - + 

s u su 



21 



= 2 -rf g 2 g 2 — [vl + 1 - 2x ffie (1 - x ee - v B )] . 



N c 



(7.117) 
(7.118) 



Counterterms 

According to eqs. (5.136), (5.145) and (5.146) of ref. [27], the counterterms are given by 

- (l + ^ee) 



1 1 



U Xa 



■2g 2 C F 



N r V ffi 



J)19,1 



1 1 



1 x ®e 
2 



\Mb(x 9B fc ,fc e )| 2 



X r 



t X 



■2g 2 C F 



SO 



C F 2 2 1 

= 2 — 5 5 S - 



1 - Xg 

2 

1 Xmc 



- (!+^ee) 



|-M,e(A; e ,x ee A: e 



(7.119) 



(7.120) 
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for gluon radiation from a qq initial-state, and 

1 1 



V 9i,g = 2g 2 s T F [l-2x 9e {l-x 9e )]\M B (x 9e k @ ,k £ 



U X a 



2 ^ g 2 9 2 s l~ [1 - 2 x ee (1 - x ee )] (7.121) 



= -IJ-2 5s 2 r F [l-2x ffie (l-x ffie )]|7W i? (A ;0 ,x a3e /c e 



= 2 -f- <? 2 9t [1 - 2 x ffie (1 - x ee )] (7-122) 
for the and (75 initiated processes, respectively. 



Radiation variables and inverse construction 

Following section 6.4.2, given the rapidity Y, the invariant squared mass of the produced 
boson M 2 and the three radiation variables (x ffie , v @ and (f>), we can construct the full 
2-body event. 

For definiteness, we assume that the emitting parton is the © one. A similar procedure 
can be followed if the emitting parton is the © one, with the exchange of the roles of the © 
and © directions. Using eq. (7.99), we can compute x ffl and x e , and from eqs. (6.91)-(6.93), 
we have 

fe ffi = x ffi K = — K, k e =x e K = x e K . (7.123) 

In the centre-of-mass system of (k B + k e ), the energy k 2 of the final-state parton and its 
angle 6' 2 with respect to the k 9 direction are given by (see eqs. (6.96) and (6.97)) 



4° = (1 " *ee) , cosfl^ 1 2 *» ^ , (7.124) 

so that 

k 2 = k' 2 ° sin 6' 2 cos <p , k' 2 2 = k 2 sin 9 2 sin <j> , k' 2 3 = k 2 cos 9' 2 . (7.125) 

We can now boost back in the laboratory frame the momentum k' 2 , with a boost parallel 
to the © direction, with velocity given by 

P = Xffi ~ Xe , (7.126) 
x ® + x e 

and obtain the radiated parton momentum fc 2 . There is no need of further Lorentz boosts 
to obtain k\ , since we can simply use momentum conservation 

k 1 = k (S + k e -k 2 . (7.127) 

According to eq. (6.90), the infinitesimal phase-space volume, d&2, can be written as 

dx 9 dx e d<5> 2 = d*i d^fad (7.128) 
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where 



d<5> 



(2A; e • k e 



rad 



16vr 2 2vr 



— dv fl 



d,X a 



@( x ®eO- -^ee)) @( x ®e x <b) • 

(7.129) 

The expression for the infinitesimal phase-space radiation volume, in the case of emission 
from the e parton, d&f ad , can be obtained with the interchange v§ <— » v e , x ffi <— » x e . The 
functions appearing in eq. (7.129) have a simple physical interpretation: they constrain 
the energy in eq. (7.124) to be positive, the cosine of the 0' 2 angle (see eq. (7.124)) to be in 
the appropriate range and force x B and x e to be in the physical range. 

Collinear remnants 

The collinear remnants of eq. (2.108) are given by 



W*i,e) = ^ 



2 (1 - zf 
log 

1 — z z 



(l + z)log^ J - + (l-z) 



+ 



+ -7T 2 - 5 ]S(1 - Z) + 



1 + z 2 

1-2 



21 



1 M 

log— n- 



B qq (zk 9 ,k e ) , (7.130) 



2vr 



T f [z 2 + (1-z) 2 ] 



, 1-z 2 M 2 

log h log — r 



+ 2z(l-z)|B, ff (zfc e ,fc e ) • 
(7.131) 



The other two collinear remnants ^e 9 (*i jQ ) and Q^{<&i^ e ) are equal to £?e 9 (3>i,e) anc l 
£?e'(*i,e) respectively, with B qq (zk 9 ,k e ) replaced by B qq (k 9 , zk e ). 

POWHEG ingredients 

We define the contributions to the real differential cross section, and corresponding coun- 
terterms, according to the three singular regions, by attaching the flux factor l/(2s) = 
x ffie /(2M 2 ), to the squared matrix elements 



w*(*i.*£d) 



9 l-M.,-1 2 , 



w«(*i.*Sd) 



2M 2 

x me 
2M 2 



|M, 9 | 2 , 



rad ~~ 2M 2 991 ' 



and 



C ffi -C*i $ ffl ,1 = V q9 ' q 

^qqV*li ^ra,&> 2 M 2 

C $ e ) — X ®° V qg ' q 



(7.132) 
(7.133) 
(7.134) 

(7.135) 
(7.136) 
(7.137) 
(7.138) 
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We also define 



99 



^qq ' qq 

C e - 
^qq 



rad 



r®- + C e - 

^qq ' qq 



(7.139) 
(7.140) 



where the counterterms are evaluated at the values specified by the radiation variables. 
The expression for B qq -(J&\) of eq. (4.13) is then given by 

=B gff (#i)+V 9ff (*i) 

+ / ^fad [*§,(*!, " ^l^fad)] 

- [G*(#i,e) + (# 1>e )] + / - [G«(*i,e) + Gf (*i,e)] , 

(7.141) 



+ 



where 



= £ !? -(S 9 ,ie) - M*0 = V«(*i) A 9 -(x e ,x e ) , (7.142) 

and the quantities R and C are obtained multiplying the corresponding quantities TZ and 
C by the corresponding factor C (x @ ,x e ). Furthermore 



<^(*i,e)=^(*i,e)4ff(*e,^), 
Gf(*i,e) =5f(*i,e)^« (*e,^f), 



(7.143) 
(7.144) 
(7.145) 
(7.146) 



with 



Cf f .(x <B ,x e )=ff(x (6 ,i4) ff,(x e ,(4)- (7.147) 
where f@{x,n%) is the parton density function of the parton / in the hadron ©, and /i 2 , is 



the factorization scale. 

Generation of the Born variable Y 

From eq. (7.129), we can write 

M 2 



rad 16vr 2 







2vr 



dx a 



dv 6 



( x ee) 



2 ' 



(7.148) 



- 73 - 



and with the change of variables 



<p = 27T(p, v @ = (1 - x ee ) v , x ee = x© + (1 - x @ )x, (7.149) 

we have 

^fad = 7^2 (1 - %)* / <W dx(l-x) dv , (7.150) 

ib7r JO JO Jo (x© + (1 - x@) x) 

so that now both d$f ad and ci$® ad are expressed in terms of the same integration variables 
and integration ranges. 

The set of variables {<p, x, v} are equivalent to the three -X" rad variables, introduced 
in section 4.3.1. We can now insert the process-dependent part of the Jacobian into the 
integrand, and define cZ 3 X rad as 

»1 rl fl 



f d 3 X iad = f dip [ dx f dv, (7.151) 
J Jo Jo Jo 



so that 



Ju[2 i 

d$? ad = d 3 X Tad J®, j*= (i-s _ (7.152) 

lovr ( x@ + (1 - x@ ) x) 



In addition, we make a change of variable in the integrals of the two collinear remnants of 
eq. (7.141), 

dz 
dx 



z = x @ + (1 — x@) x, and — = 1 — x@ , (7.153) 



so that the limits of x are consistently between and 1. The B function of eq. (4.23) is 
then given by 

5 M -(*i,X rad ) = + 

+ J ffi [i? fl? -(* 1 , $? ad ) - C gq (^, + J e [R qg (^ - $ r e ad )] 

+ ^ [G?(*i,e) + GT(*i,e)] + ^ [GT(*i,e) + G?(*i,e)] , 

(7.154) 

so that 

J d 3 X rad B M -(*i, X rad ) = B M -(*i) . (7.155) 

The generation of the Born variable Y is performed using the method illustrated in sec- 
tion 4.3.1. We use an integrator-unweighter program (like for example the BASES-SPRING 
package), that, after a single 4-dimensional integration of the function 

B(*l,X rad ) = 2^ M -(*l.*rad) , (7.156) 
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can generate 4-tuples of {Y, ip, x, values, distributed according to X rac j) . For 

each generated 4-tuple, we generate q with a probability proportional to B q q(&i, -X" ra d) - 19 
We only keep the Y and q generated value, and neglect 92, x, v, which corresponds to 
integrate over them. 

Observe that B(&i, X rax \j is positive, unless the 0(a s ) terms overcome the Born term. 
If this happens, the whole perturbative approach breaks down. This can happen if M is 
too small, or if we are in extreme regions of phase space, like the threshold region (i.e. 
when the values of Y and M are forcing x ffi or x e to 1). In practical applications, it is wise 
to check that the fraction of negative weights in the integral of B(&i,X ra d) is small, and 
can be safely neglected. 

Generation of the radiation variables 

The transverse momentum k T of the radiation, with respect to the incoming beams, in the 
two singular regions is given by 

k 2 T@ = s(l-x ee )d® = M 2 ^^v®, (7.157) 

x s>e 

where is used when the emitter is the e parton, and v e when it is the e one, and the 
Sudakov form factor of eq. (4.16) is 



-/ 



In order to generate the radiation variables, we use the hit-and-miss technique introduced 
in section 4.3.2. Suitable upper bounding functions are given by 

*^rad) <- o (~i „2 X ffi9 J_ ^ /<£ { x ®/ x ®e ; Mf) /y i rn\ 

-^<?q(^li ^rad) x fi( o „2 X S9 ^ ^ fe ( X e/ X <BQ i Mf) / 7 i Kn N 

ggg(£ii^|d) < OT 2 ^ 1 fj(xjx @e ,i4) , 71fin 

Rqg(&l, ^fad) <r o T „2 x ee 1 /e(^e/ x ®9) Mf) / 71fio \ 

^(♦o - 8Tf9 '2^«; ici^i) ' (7 6) 

Following appendix D, we can put an upper bound on the ratios of the parton distribution 



The standard procedure to generate an index < j ^ n with a probability proportional to a,j is to 

^"Jj ay < r £™, =1 ay < Ey=i ay. 



generate a random number < r < 1, and then choose j so that Yjj'=i a j' < r Ey=i a j' < Sy=i a 
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functions (pdf) too, and we obtain 



y y v 1 ' raa / 
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(7.163) 
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(7.164) 




1 






^-(*l^fad) 


< AT-^ 
u ® 




X 


Be 


(7.165) 
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^©e 
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Be 


(7.166) 




1 




^©e 



where N^, N g q and N qg are upper-bound constants determined by sampling the whole 
phase-space of the R/B ratios. 

We can then write an upper bound for the integrand function in the Sudakov form 
factor of eq. (7.158) 



/ 



d<5> 



E>© I O 



J Dqq 



rad 



B. 



qq 



< iV / dx 0e / ete® tt^ r — 

Jx® Jo x ©e I 1 x ©eJ v © 

+ iV e / cfa ee / 

Jx JO 



«s (fe? e ) ^(^t® -Pt) 
a s (fe? e ) 0(*£ e -p£), (7-167) 



/z e jo x @e (1 x ©e) 

where N & are constant with respect to the radiation variables. The two integrals have 
the same functional form, so that we just need to generate radiation variables using the 
highest-p T bid method of appendix B and the veto technique of appendix A, where the p T 
is generated uniformly in A u (p T ) 



Au(Pt) = exp 



-N 



[ dx [ 

Jx JO 



rir — t~~ r- a s (k 2 ) 6(kl - p%) \, (7.168) 



x (1 — x) V 

where k T must be of the order of the radiation transverse momentum in the collinear limit, 
and it must coincide with it in the soft-collinear limit. A suitable definition is given by 

k% = M 2 {l-x)v. (7.169) 

Introducing the dimensionless quantities 



krw 



J; = (l-*)». 



p% 

M 2 ' 



(7.170) 



and the integration over dk 2 , we have 



A u (Pt) = exp • 



exp • 



exp • 



-N f T ' maX d^ f dx f 
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J p\ J X 
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ft rp 



log 
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+ log 
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(7.171) 
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where 



(7.172) 



and a s (fc 2 ) is given by eq. (7.45). Since the generation of p T according to 



dA u (p T ) = dexp < — N 



fc T,ma X , 2 a s (& 2 ) 



dki 



&t 



1 — k T 1 — x 

log — = h log 



0(kl- P 2 T ) 



(7.173) 

is still too complex to be performed analytically, we use a second time the veto method, 
using the upper bounding function 



i 1 !- fc T 
log -=- > log — ~ — , 

and we then generate the p T distribution according to dA uu (p T ), where 



(7.174) 



A uu (Pt) = exp 



cxp 



iV 
+ log 



~2 *-^s (^t) 
1,2 

fbrp 

1, M 2 (1 - x) 2 



log J- + log 1 



X 



X 



9{K-pl) 



log((l -x) 2 M 2 /A 2 ) 
log (p 2 /A 2 ) 



, 1 - x 1, M 2 
log — + - log —j- 

x 2 A A 



(7.175) 



In order to apply the highest-p T bid technique to eq. (7.158), we need to repeat the following 
steps for both the last two integrals in eq. (7.167). We follow in detail the first one. The 
second integral is treated in the same way. 

Observe that we have a lower limit on the acceptable values of p T , since a s {p T ) must 
be well defined for both the two- loop and the one-loop expression of a s . We thus introduce 
a p™« > A. 

When dealing with the first integral, we set N = N 9 and x = x ffi . 

1. Set p™ ax = M (1 — x), where M 2 (1 — x) 2 is the maximum value of p 2 , such that 

2. Generate r with < r < 1, and solve the equation r = A uu (p T ) / A uu (p™ x ) for p T . 

If no solution is found, or if p T < p™ n , no radiation is produced, and a Born event is 
generated. 

3. If a solution with p T > p™ m is found, generate r' with ^ r' ^ log(M/p T ). If 
r' ^ log(M/p T — 1) go to step 4. Otherwise set p™ ax = p T and go to step 2. 

4. At this point p T is generated according to eq. (7.173). We need to generate x, v and 
(f) with a probability proportional to 



1 



1 



x (1 — x) V 



-5{M 2 v{l - x) -p 2 T )dvdxd(t). 



(7.176) 
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This is done as follows: integrate the above in v using the 6 function, so that 



vl 1 



(7.177) 



M 2 1 - x ' 
and generate x according to 

e ix -x) e(l-^-x) dlogr^, (7.178) 

i.e., generate a uniformly-distributed random number r" with 

log < r" ^ log (— - l) (7.179) 
1 - x \p T J 



and solve r" = log (x/(l — x)) for x, that is 

1 



exp(— r") + 1 ' 

One also generates uniformly a random value of 4> between and 2ir. 



(7.180) 



5. Set x ee = x and v B = v. Now we need to apply the last veto. Generate a random 
number r ffl in the range 

0<re<W + Ay^T^. (7.181) 

^(•..O + JW*..*^ (7 . 182) 

then accept the event. Otherwise, set p™ ax = p T and go to step 2. 

6. Call the p T generated this way p®. Set p® = in case the event is returned as a 
Born one. Repeat the previous steps for the second integral, that is set N = N e and 
x = x e , with the obvious changes in step 5, that now becomes 



5'. Set 

x <se — x an d v e — v and generate a random number r e in the range 

If 



0<r e <(N^ + N qg )^-^. (7.183) 



then accept the event. Otherwise, set p™ ax = p T and go to step 2. Call the p T 
generated this way p®. Set p®. = in case the event is returned as a Born one. 

7. According to the highest-p T bid method of appendix B, choose the highest value 
between p® and p®, . If both p® and p®< are zero, then return the event as a Born-like 
one. 
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• In case the p T chosen is p® , if 

then generate a gq event. Otherwise, generate a qq event with radiation emitted 
by the e incoming quark. 

• In case the p T chosen is p® j if 

then generate a qg event. Otherwise, generate a qq event with radiation emitted 
by the e incoming quark. 



7.4 h e h G — > V in the Frixione, Kunszt and Signer formalism 

We now consider the case of the production of a massive vector V of mass M in hadron- 
hadron collisions, in the FKS formalism. 



Radiation variables and inverse construction 

We use the kinematics notation introduced in section 7.3, that is, the Born kinematics is 
characterized by the rapidity Y of the produced vector boson, and the phase space and 
momentum fractions are given by eqs. (7.98) and (7.99) 




(7.187) 



According to eqs. (7.105)-(7.107), the real-emission processes are characterized by two 
final-state momenta, k\ (the V momentum) and k 2 , the momentum of the radiated parton. 
The (only) radiated parton is the FKS parton. Following the prescriptions of section 5.1, 
we define 

fs 

fc° = ^-£, k 2 = ££(l,sin6>sin^,sin6>cos<£,cos6>), (7.188) 

and we take £, y = cos 6 and 4> as the radiation variables. The real-emission phase space is 
given by (see eqs. (5.9) and (5.10)) 

d# 2 = d*i d$ rad , d$rad = t/to zr^-z <% dy d<P , (7.189) 

(47T)° 1 - £ 

and the procedure described in section 5.1.2 allows one to fully reconstruct the two-body 
phase space, give the Born barred variable Y and the three radiation variables. In partic- 
ular, the requirement that both x ffi and x e of eq. (5.11) be less than 1 forces £ to be in the 
range < £ < £ max , where £ max is given by eq. (5.13). 
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Virtual corrections 

The virtual corrections to the process can be found in ref. [37], eq. (82) 



Vb,qq = 



47r/x 2 y r(i + e )r(i - e ) 2 « s 

M 2 ) r(l - 2e) 27 



C v 



2 3 

8 + ^ 2 

e 



(7.190) 



where we drop the argument <&i of the virtual and Born term, for ease of notation. Fol- 
lowing section 7.3, we indicate with q the Born flavour index fa, that carries information 
on the flavour of the (massless) incoming parton, along the © direction. Using a standard 
convention for the numbering scheme, we have q = —5, . . . , —1, 1, . . . , 5, and q = —q. We 
notice that 

r(l + e)r(l-e) 2 _ 1 



r(l - 2e) 



r(i-e) 



+ O (6 3 ) , 



so that, choosing Q = M in eq. (2.93), we get 



Vb,qq 

From eq. (2.92) we thus get 



e 2 e 



B 



qq- 



V fin , w - = C P (7r 2 -8) B, 



qq ■ 



(7.191) 



(7.192) 



(7.193) 



Notice that the Bij terms in (2.92) have only the £> ee and £> e0 components in this case, 
and their contribution vanishes because we have chosen Q 2 = M 2 = 2 (k 9 ■ k e ). According 
to eqs. (2.97), (2.100) and (2.101) we have 



B a 



B, 



C^Bqq , 



(7.194) 



lij Bjj = (B^q +B e(B )l me = -(B Be + Z3 effi )Li 2 (l) 



TT 



Q = -2C F log 



/4 
M 2 ' 



— (B ee + H ee ), (7-195) 



(7.196) 



where we have chosen £ c = 1. Finally, eq. (2.99) gives 



qq ■ 



(7.197) 



Collinear remnants 

From eq. (2.102) we immediately get the collinear remnants 



as 



1 



1 - z 



, M 2 
log — + 2 



log(l - z) 



1 - z 



+ 1 Zf Bqq(zk <£ , k e ) , 

(7.198) 



Qf=^T F \[z^ + {l-zf] 



M 2 

log — + 21og(l-z) 

zn 2 



+ 2z(l-z)^Bqq(zk 9 ,k e ),(7.199) 
and the corresponding ones for the e collinear direction, and where we have set S t = 2. 
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Real corrections 

The real production contributions are given by eqs. (7.132)-(7.134) (see also [37]) 



1 

2~s 



K 



fi'i 



r> 2 2 



N, 



2s 



u t 2 s M 2 

t u tu 

u s 2tM 
- + - + 

s u su 



21 



(7.200) 
(7.201) 



(and the analogous one for 1Z qg ) where, according to the definitions given in eqs. (7.108)- 
(7.110), 



M 2 s _ . s , . 

s= r3^' t = --e(l + y), u = --C(l-y) 



(7.202) 



In terms of the FKS variables we have 



(i-0 



Cp 2 2 1 [ l+j/ 2 , 

'^- 2 N C 9 9s 2s\ 1-y 2 + eil-y 2 )^ 
v T F 2 2 1 f 2[l-g(l-fl(l + y)] g(l-y) 



2s 



(7.203) 
(7.204) 



so that (1 — y )£ H are regular functions. Using eq. (2.84) we obtain immediately 



-4 ^ J. [2 (i + n e + 8( i - {)] {i (i) [(^ 



+ 



^- = 2^5 2 5 S 2 ^^2[1-<(J-0(1 +//)] + 



e(i-y) 2 \ ( i 



l 

!' 

(7.205) 
(7.206) 



with the + distributions defined in eq. (2.90). 

Generation of the Born variable Y 

The expression for B q q(&i) of eq. (4.13) is then given by 



B M -(* 1 ) = 5 M -(* 1 )+^-(*i) 



+ J d$ lad flg ? (*i,$rad) +-Rg fl (*l,^rad) + -Rgf (*1 ' $ rad) 

f 1 rlz - - - - f 1 Hz 

+ I - [Gg*(#i, e ) + GT(*i,e)] + / - [G?(*i,e) + Gf(*i,e)] , (7-207) 
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where 



i?qq(*l,$ rad 

G|f(*i,e 
Gr(*i,e 



#gg(*l) £gg (x ffi ,X e ) , 
Vgg(*l) £gg (x e ,X e ) , 

K q q(&i, $ rad ) £ 9 g (ar e , x e ) , 
n qg (&i, ^rad) £ 99 ( x e, a?e) , 

tigqi&l, $ ra d) (^e.^e) > 

Sf(*i > e)A ff f^,x e ), 



z / 



^(*i,e)^ 
af(* 1)e )£ ?9 -(x e ,^), 

0I ff (*i,e)£ w 



(7.208) 
(7.209) 
(7.210) 
(7.211) 
(7.212) 

(7.213) 
(7.214) 
(7.215) 
(7.216) 



with x ,x e given in eq. (5.11) and the luminosity C defined in eq. (7.147). 

As explained in section 5.1.1, the integration range in the radiation variable £ is limited 
by the condition (5.12) 

< e < 6nax • (7.217) 

Making the change of variable 

l-£ = z, (7.218) 

the integration limits become 

Zmin(y) < Z < 1 , (7.219) 



where 



^min(y) = max 



2(1 + y)* 2 



V(i + 4) 2 (i-y) 2 + i6y4 + (i - 
2(1 -y) x| 



37 /T 



V(i + 4Hl + yp^l6y4 + (1 + y)(l - x%) 
It is convenient to introduce a rescaled variable z 

Z ^-min(y) 



^min (2/) + « (1 - «min(y)) , 



(7.220) 



(7.221) 



and use the identity 



1 



1 



1 



{ ( rb) + 5(1 _ 5) log[1 " Zmin(y)] } ' (7 - 222) 



in order to simplify the treatment of the + distribution in the integration of R. 

The integral of the distributions requires some care. We follow, as an example, the 
integral of the 1/(1 — y) + contribution in the R qq - term. The integral has the form 



f dy f 

J-i Ji-, 



1 

i^y 



+ 



f(Z,v)- 



(7.223) 
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We change variable from £ to z using eq. (7.222), and obtain 

+ J_ 1 dy (rh;) log [1 " '™ |,<!/)I f{0 - y) 

log [1 - z miQ (y)] f(0,y) - log [1 - z min (l)] /(0, 1)} , (7.224) 

that is manifestly finite and can be computed numerically. 

The collinear contributions also have a restricted range in the z variable. Considering 
for simplicity only the term, we have x B < z < 1, so that, as usual, we need to make 
use of the identities 

j\z(-?-) f(z) = log(l - x»)/(l) + jf ' <fc /( '| ~ ^ W , (7.225) 

/' J(Ml^) + = 'vp -,.,/d) 7/' ^MMiM, ( , 226) 

Ji® V 1 Z /+ Z J% 1 Z 

and also in this case we define a rescaled variable < z < 1 

. = Z __I®. ; z = 2e> + 5 (i _ j ( 7 .227) 

i x® 

and rewrite the z integrations as z integrations. At the end of this procedure, the most 
general form one can obtain for B is 

rl rl r2n 

B M -(*x) = B a qq (^) + / d5B;,-(*!,5) + dy d<j>B c qq (* u yA) 

JO J-l JO 

+ / dy / dz / d4B* q {* u z,y,<t>) . (7.228) 
7-i Jo jo 

In order to perform numerically the generation of the Born term, we introduce the function 
(see eq. (4.23)) 

B qq (* 1} z,y,<i>) = ^B^ 1 ) + ^-£dzB b qq (^ 1 ,z) + j\y J** fyB^^y^) 

/l rl r2n 

d vJ d *J <tyB*-(*i,z,y,<£), (7-229) 

so that 

/l rl rln 

dy dz d<t>B qq {^ u z,yA) ■ (7.230) 
-1 jo jo 

The generation of the Born variable <&i = {Y} is performed by using an integrator- 
unweighter program (like for example the BASES-SPRING program), that after a single 
four-dimensional integration of the function 

z,y,<f>) z,y,0) (7.231) 
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in the variables Y, z, y, 4>, can generate 4-tuples of Y, z, y, (f> values, distributed according to 
the weight B(Y, z, y, (f) . For each generated 4-tuple, we also generate q with a probability 
proportional to B q q(Y, z,y,(f)). We only keep the Y and q generated value, and neglect 
y, z, cf>, which corresponds to integrating over them. 

As already noticed at the end of section 7.3, observe that B(Y,z,y,4>) is positive, 
unless the 0(a s ) terms overcome the Born term. If this happens, the whole perturbative 
approach breaks down. This can happen if M is too small, or if we are in extreme regions 
of phase space, like the threshold region (i.e. when the values of <&i and M are forcing x 9 
or x e to 1). In practical applications it is wise to check that the fraction of negative weight 
in the integral of B(Y, z, y, cf) is small, and can be safely neglected. 

Generation of the radiation variables 

We assume that we have generated <&i and q according to the procedure given above. The 
Sudakov form factor for the generation of radiation, according to eq. (4.16) is given by 

A /=== v f /",_ flrc(*l,*rad) +fl M (*l,$rad) 1, $rad) an \ 

Ag(*i >Pr ) = exp<-y d$ rad ^ ^ '- 9(k T -p T )>, 

(7.232) 

where 

kl = {e{l-y 2 ) = ^ l) e{l-V 2 ) (7.233) 

is the exact squared transverse momentum of the radiated parton. The factorization and 
renormalization scales in eq. (7.232) should be taken equal to k%. In order to generate the 
radiation variables, we use the hit-and-miss technique introduced in section 4.3.2. We need 
to find a simple upper bound for the expression 

8 £ Rqq{*l, $rad) + Rqg{*l, $rad) + Rgq(&l, ^rad) (? ^ 



(4tt)3 1-£ S m -(*i) 

in order to perform the generation of radiation using the veto method. With the same 
considerations of section 4.3.2 and using the upper bounds for the parton distribution 
functions found in appendix D, it is easy to show that the upper bounding function found 
in ref. [22] 

U q = N q ^0L (7.235) 

is suited also to the present case. 

When implementing the bound (7.235) in practice, the whole phase space is searched 
randomly in order to determine the bound normalization iV^. At this stage, problems 
can arise, due to the fact that available parton-density parametrization do not necessarily 
comply with the bound (D.8). One simple example is the case of the b quark density, that 
is set arbitrarily to zero above a given value of x, generally smaller than the value of x at 
which the gluon density is set to zero. These regions give very small contribution to the 
cross section, and must be carefully cut away, unless one is using a pdf parametrization 
that is fully consistent with evolution also in the very large-x region. 
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The bound (7.235) has been obtained assuming that we are not near the small- x region. 
If this is the case, the bound may become inefficient. We observe that, in case the parton 
distribution functions obey Feynman scaling (i.e. f(x) « 1/x at small x), the bound is still 
adequate in the small-x region. Modern pdf parametrizations, however, prefer the faster 
rising behavior f(x) « 1/x 5 , with 5 > 1. In this case, the bound may become inefficient 
at small x, and numerical tricks (like, for example, dividing the integration range in small 
intervals and adopting different N q values in each interval) should be used to overcome this 
problem. 

The generation of the event according to the bound (7.235) is documented in great 
detail in appendix D of ref. [22], and we do not repeat it here. 

8. Conclusions 

In this work we have presented in full details the POWHEG method for interfacing NLO 
calculations to parton shower. The purpose of this paper, besides providing all the technical 
details needed for the implementation of the method, was also to demonstrate its full 
generality. We have given detailed formulae in the case of cross sections that do not 
involve massive coloured partons. The extension to this case is straightforward both in the 
FKS and in the CS framework. We have not tried to include it in our discussion, in order 
not to increase the complexity of the presentation. 

We believe that the method presented here can be applied to a large number of inter- 
esting LHC processes, and that authors of QCD calculation should be able to build their 
own implementation with only a modest effort. In order to ease this task, we have collected 
POWHEG software in a repository at the location: 



Among other things, the software for the implementation of the inverse mapping in FKS, 
the program mint (see ref. [38]) for the integration and generation of unweighted events, 
and the full implementation of the e + e~ example in the FKS approach can be found there. 
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A. The veto technique 

In this section we describe a method to generate a set of d variables x, distributed according 
to 



http : //moby .mib. inf n. it/~nason/PDWHEG/FNDpaper/ . 



f{x) A(h(x)) 



(A.l) 



where 




(A.2) 
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We assume, as usual, that / and h are non-negative functions, and that the unrestricted 
integral of / is divergent, that is 

A(0) = expl- J d*a/f(a/)\ =0. (A.3) 

With these assumptions, upon multiplying the infinitesimal probability f(x) A(h(x)) d d x 
by 5(h — h(x))dh, we can integrate over d d x and interpret it as infinitesimal probability 
for the variable h 

dh J d d x5(h-h(x))f(x)expl- J d d x'f(x') 9(h(x') - h)\ = dh^^- = dA(h) , 

(A.4) 

that shows that the probability is uniform in A(h). In principle, the generation of events 
is therefore straightforward: one generates a uniform random number r between and 1, 
and solves the equation A(H) = r for H (here we have used the fact that A(0) = 0). At 
this point, the variables x have distribution function equal to 

6(H - /i(x))/(x)exp|- J d d x'f(x') 0(h(x') - H) j , (A.5) 

where the exponent is now just a number (a normalization factor), so that the variables x 
are on the surface 5{H —h{x)) with a distribution function proportional to f(x) 8(H—h(x)). 
The generation of these variables can be done, for example, with a hit-and-miss technique, 
or, if the integration can be performed analytically, by generating (d— 1) random numbers 
rj, uniformly distributed between and 1, and solving 

[ X * dxi H f dx k 5{H - h(x)) f(x) = n (A.6) 



for Xi, where Xio is the lower limit for the variable Xi, and all the other variables are 
integrated over their full range of validity. 

In practise, however, the solution of the equation A(H) = r is, in most cases, very 
heavy, from a numerical point of view. This difficulty can be overcome by means of the 
so-called veto method. We assume that there is a function F(x) > f{x) for all x values, 
and that 



Ap(h) = exp 



d d x'F(x') 9(h(x') -h) 



(A.7) 



has a simple form, so that the solution of the equation Ap(H) = r and the generation of 
the distribution F(x) 5{H — h(x)) are reasonably simple. Then, we implement the following 
procedure: 

1. Set H mSLX equal to the maximum allowed value, such that AF(H mSLX ) = 1. 

2. Generate a flat random number < r < 1 and solve the equation 

A H#) = r (A8) 

Ai^i^max) 

for H (a solution with < H < H nmx always exists for < r < 1). 
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3. Generate x according to F(x) 5(h(x) — H). 

4. Generate a new random number r' . 

5. If r' > f(x)/F(x) then the event is vetoed, we set H max = H, go to step 2 and 
continue. Otherwise the event is accepted, and the procedure stops. 

The resulting events are distributed according to eq. (A.l). The proof of this statement is 
simple but non-trivial, and is given in appendix C of ref. [22]. 

B. The highest-p T bid procedure 

Our aim is to generate (k, Xk) pairs with a probability 

f k (x k )l[\(h k (x k )) d d x k , (B.l) 

i 

where 

A i (h) = exp{-J d d x' % U{x[) B{hi{^-h)Y (B.2) 

We assume, as usual, that the fi and hi are non-negative functions, and that the unre- 
stricted integral of the fi is divergent, that is 

Ai(0) =exp|- J rf^/i^) j =0. (B.3) 

Under these conditions, we have the identity 

/d 

d d x k A k (h k (x k )) fk(x k ) S(h k (x k ) - h) = ^ A k(h) , (B.4) 

so that 

/A./Mfc(*»))A(».) «(*-»»(»»)) 

poo r 

= dh' J d d x k A k (hk(x k )) fk(xk) 6(h' - h k (x k )) 0{h- h') 

poo j fh j 

= ^ dh' -^Ak{h') 0{h - h') = dh' —A k (h') = A k (h) , (B.5) 

where we have used the fact that A&(0) = 0. If we interpret h and hk as transverse 
momenta, then Ak{h) in eq. (B.5) corresponds to the probability of not emitting a parton 
with transverse momentum bigger than h. 

The procedure to generate the distribution in eq. (B.l), using the highest-bid method, 
is the following. For each k, we generate an Xk value with probability 

Ak{hk{x k )) fk(xk) d d x k , (B.6) 

as described in appendix A, and then we pick the k value with the largest hk(xk)- In fact, 
the probability that the generated (k, x k ) has the largest hk{xk) is precisely given by the 
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product of its generation probability (eq. (B.6)) times the probability that all the other 
hi(xi) are less then hk(xk), which is given by the product 

Y[Ai{h k (x k )), (B.7) 
and together they reconstruct eq. (B.l). 
C. Soft angular integrals 

The soft angular integrals eqs. (4.45) and (4.46) can all be obtained from the basic integral 

'^'•^ = hwm^ry (ai) 

where v is an arbitrary timelike vector and dfl is the element of the solid angle of k with 
respect to a reference direction. First we notice that, since the integrand scales like (k°)~ 2 , 

f r e f( n ) h. . „ 

dn k°dk° 1 =i(ki,v,i), (c.2) 

7 7/(n) (ki-k)(k-v) 

(where e is the Euler constant) independent of /, for any arbitrary function /(fi). We can 
also write 

m,v,D = j ^ {ki k k ){ v k . v) Hk°-m)) e(em)-k°), ( c. 3) 

and perform a change of variables, that is also a Lorentz transformation, k — > Afc, to obtain 

Kh, v,l)=J ^ ^i^ - /'(")) '(*/'(«) - kP) , (C.4) 

where v' = A~ 1 v and k\ = Ar l ki. The function / also undergoes a change, and becomes 
/'. But we have just shown that the integral does not depend upon /, so that we conclude 

i(ki, v ,i) =/(^y,i). (c.5) 

We thus perform the integral in a frame where v has only the time component. We get 



Setting v = h L j in eq. (C.6) yields immediately eq. (4.46). For small kf, eq. (C.6) becomes 

tii iO\ 27r i ^{h-vf 

I(ki,v,k u ) = -^\og fc2 ^ 2 . (C.7) 
Taking the difference of the above formula with v = ky and v = n we get (4.45). 
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D. Parton-distribution-function upper bounds 



In this section, we discuss the bounds on the parton-distribution-function (pdf) ratios of 
the form 

fg(x/z,/4) f q (x/z,l4) f g (x/z,l4) fq(x/z,f4) 



fg(.B,l4) fq(x,(4) ' fq(x,(4) ' fg{x,f4) 



(D.l) 



where q stands for a quark or an antiquark, and x < z < 1. We begin by noticing that the 
first two ratios are always bounded by a number of order 1, because the pdf's are never 
fast growing functions of their argument. In fact, only the pdf of the valence quarks grows 
mildly for moderate x values. In the remaining two cases, the ratio of parton densities is 
between different parton species, and must be discussed with care. We consider first the 
g/q case. First of all, since the gluon pdf is a monotonic decreasing function 

f g (x/z,i4)<fg(x,i4), (D.2) 

so that we only need a bound for 

fg(x, fi F ) 

At small values of x, this ratio is always bounded, because the gluon and quark densities 
have similar behavior in the small- x limit. For large x, if q is a valence quark, the ratio is 
also bounded, since the gluon is generated by valence quarks through evolution. In case 
quark, the corresponding density may be softer than the gluon density. In the 
worst case, however, the sea quark is generated by the gluon through evolution. Assuming 
that parton densities behave as a power of 1 — x at large x, 

f g (x,£)~(l-x) s , (D.4) 

the Altarelli-Parisi equation in the large-x limit yields 

2 ^4) ^ TVas /* dz ^ JVas s+l 



and therefore 
Since 

we conclude that 

For the q/g case, since 
we need a bound for 



fifed® < -°- . (D.6) 

fq(x,l4) ~ 1 - X 

1 <T 1 T> (D.7) 



1 — x 1 — z 

fg(x/z,(4) C 

-jj^w - — ■ (D - 8) 

fq(x/z,l4) < fq{x,(4), (D.9) 



2^ 



fq(x,H 



(D.10) 
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The symbol < in eq. (D.9) is to be interpreted as "almost everywhere smaller than". In 
fact, if q is a valence quark, its pdf can grow mildly in an intermediate x range. If q is not 
a valence quark, the ratio (D.10) is certainly small. We should only worry about the case 
when q is a valence quark. In this case, we assume for the quark a large x behavior of the 
form 

f q (x,£)~(l-x) 5 , (D.ll) 
and, with a reasoning similar to the one used before, we can conclude that 

f q (x/z,(4) C 
f B M) -1-z- [D - 12) 
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